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Barack and Sago [Phys. Rev. Lett., 102, 191101 (2009)] have recently computed the shift of the 
innermost stable circular orbit (ISCO) of the Schwarzschild spacetime due to the conservative self- 
force that arises from the finite-mass of an orbiting test-particle. This calculation of the ISCO shift 
is one of the first concrete results of the self-force program, and provides an exact (fully relativistic) 
point of comparison with approximate post-Newtonian (PN) computations of the ISCO. Here this 
exact ISCO shift is compared with nearly all known PN-based methods. These include both "nonre- 
summed" and "resummed" approaches (the latter reproduce the test-particle limit by construction) . 
The best agreement with the exact (Barack-Sago) result is found when the pseudo-4PN coefficient of 
the effective-one-body (EOB) metric is fit to numerical relativity simulations. However, if one con- 
siders uncalibrated methods based only on the currently known 3PN-order conservative dynamics, 
the best agreement is found from the gauge-invariant ISCO condition of Blanchet and Iyer [Classical 
Quantum Gravity 20, 755 (2003)], which relies only on the (nonresummed) 3PN equations of motion. 
This method reproduces the exact test-particle limit without any resummation. A comparison of 
PN methods with the ISCO in the equal-mass case (computed via sequences of numerical relativity 
initial-data sets) is also performed. Here a (different) nonresummed method also performs very well 
(as was previously shown). These results suggest that the EOB approach — while exactly incorpo- 
rating the conservative test-particle dynamics and having several other important advantages — does 
not (in the absence of calibration) incorporate conservative self-force effects more accurately than 
standard PN methods. I also consider how the conservative self- force ISCO shift, combined in some 
cases with numerical relativity computations of the ISCO, can be used to constrain our knowledge 
of (1) the EOB effective metric, (2) phenomenological inspiral-merger-ringdown templates, and (3) 
4PN- and 5PN-order terms in the PN orbital energy. These constraints could help in constructing 
better gravitational-wave templates. Lastly, I suggest a new method to calibrate unknown PN terms 
in inspiral templates using numerical-relativity calculations. 

PACS numbers: 04.25.Nx, 04.25.-g, 04.25.D-, 04.30.-w 



I. INTRODUCTION AND MOTIVATION 

The primary purpose of this study is to compare re- 
cent gravitational self-force (GSF) calculations of the in- 
nermost stable circular orbit (ISCO) [l|-|3( with nearly 
all post-Newtonian (PN) and effective-one-body (EOB) 
methods. The first half of this paper provides introduc- 
tory material, reviews previous related work, and sum- 
marizes the various PN/EOB approaches. Readers wish- 
ing to skip this material can proceed directly to the re- 
sults in Sec. [TV] 

A. Regimes of the relativistic two-body problem 

One of the goals of this study is to provide insight 
on the various methods used to solve the relativistic two- 
body problem for the purpose of generating gravitational- 
wave (GW) templates. We begin by briefly reviewing 
these methods. 



The post-Newtonian (PN) approximation iteratively 
solves Einstein's equations using the approximation that 
a binary's relative orbital speed v is small compared 
to the speed of light c. The PN equations of motion 
are known completer^ to 3.5PN order [i.e., computed 
to order (t> 2 /c 2 ) 3 ' 5 beyond the Newtonian terms; see @ 
for references and a review]. For a binary with masses 
mi < m-2 and v/c <C 1, the PN approach is valid for 
any mass ratio q = mi/m.2 < 1, although it is known to 
"converge" more slowly if q <C 1 043- 

When the binary separation is small and v/c ~ 1, 
the PN approximation breaks down, and other methods 
must be applied. One such method is numerical relativ- 
ity (NR), the numerical solution of Einstein's equation 
without approximation. This approach has had much re- 
cent success (see [l3l - [T6j for reviews), but computational 
limitations currently restrict it to modeling binaries with 
mass ratios q > 0.1 [17] (however, see Refs. [l8l - [20T | for re- 
cent progress). For smaller mass ratios the time to inspi- 
ral increases like Xi nsp ~ 1/g, and multiple spatial scales 
(Ar ~ W2 and ~ mi = qm^) must be resolved accurately. 
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Partial results at higher orders include the 4PN tail contribution 
U and the 4.5PN radiation-reaction terms 
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This requires a finer spatial grid, smaller step-sizes, and 
longer evolution times. It will therefore be very difficult 
for NR to simulate more than a few orbits for binaries 
with very small mass ratios (q < 10~ 2 ). 

Because they will execute many observable orbital cy- 
cles in the highly relativistic (v ~ c) regime, an accu- 
rate description of extreme (q < 10 ) and intermediate 
(10 -4 < q < 10~ 2 ) mass ratio binaries are amenable to 
neither PN nor NR methods. But they are amenable 
to a third method — the gravitational self-force approach. 
This is based on computing how a point-particle with 
mass mi -C m 2 deviates from geodesic motion around a 
black hole (BH) with mass mo. The force that causes 
this deviation (the GSF; see |2ll— [24jj for reviews and ref- 
erences) arises from the particle's own gravitational field. 
The GSF is responsible for dissipative effects like the 
radiation-reaction force that causes the point-particle to 
lose energy and angular momentum to GWs as it inspi- 
rals. It is also responsible for conservative effects which 
are time-symmetric and preserve the orbit-averaged con- 
stants of the motion. 

One example of a conservative GSF effect is the shift in 
the periastron advance angle per orbit Aip due to finite- 
mass ratio corrections: e.g., at leading PN order the pe- 
riastron advance can be written as 1251 



Aip 
~2T 



= k = 



3(m 2 n) 2/3 
(1-e 2 ) 



1 
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0(q 2 ) 



(1.1) 



where the mean motion is n = 27r/P or b, P Q rb is the 
periastron-to-periastron period, and et is the "time" ec- 
centricity appearing in the quasi-Keplerian formalism of 
[25| . The first term represents the geodesic contribution; 
the O(q) term represents the first-order conservative GSF 
correction. Conservative GSF calculations of the perias- 
tron shift are discussed in [26|, [2lJ ■ Another example of a 
conservative GSF effect is the finite test-mass shift in the 
frequency of the ISCO (which is the focus of this study). 

While evaluating the full GSF has proven to be tech- 
nically difficult, in the past three years four independent 
groups @, l28l - [32j have succeeded in computing the GSF 
for circular geodesies in SchwarzschildQ More recently, 
Barack and Sago (BS) have computed the GSF for ec- 
centric (bound) geodesies in Schwarzschild [l], [|| . They 
were then able to compute the change in the ISCO radius 
and angular frequency due to the conservative-piece of 
the GSF. This represents a significant milestone in (typi- 
cally gauge-dependent) GSF calculations since the ISCO 
shift is a well-defined and easily understood strong-field 
quantity that can be compared with other approaches. 
As previous studies [35M55| have investigated the agree- 
ment between NR and PN-based waveforms in the q ~ 1 



2 The GSF research program is strongly motivated by the need 
to produce accurate waveforms for extreme-mass-ratio inspirals 
(EMRIs), an important source for the LISA mission [33| consist- 
ing of a stellar-mass compact object inspiraling into a massive 
BH [H. 



regime, the objective of this study is to further compare 
PN-based approaches with these new GSF results (which 
are exact in the q <C 1 limit). 

While we investigate several PN-based approaches be- 
low, let us briefly highlight the effective-one- body mm 
approach, which has especially motivated this study. The 
EOB formalism attempts to improve the convergence of 
the PN two-body equations of motion by mapping the 
PN two-body Hamiltonian for the masses mi and m 2 
to an "effective" Hamiltonian that (at the 2PN level) de- 
scribes a particle with reduced mass fi — mim 2 /M = rjM 
moving on geodesies of a "deformed" Schwarzschild back- 
ground associated with a mass M = mi + m 2 (At 
the 3PN level the effective Hamiltonian must be sup- 
plemented with additional terms that do not arise from 
the Hamiltonian of the effective metric [60].) The ef- 
fective Hamiltonian describes the full conservative dy- 
namics of the (/x, M) binary, while quantities like the "77- 
deformed" ISCO, light-ring, and horizon depend only on 
the time-time piece of the EOB effective metric function 
g^f = —A(r). To include the effects of dissipation, the 
EOB approach incorporates information from the PN ex- 
pansion of the energy flux (resummed by various means 
to improve convergence). Using these elements, the EOB 
formalism is able to describe not only the inspiral, but 
also the transition region where the inspiral ceases to be 
adiabatic and the point mass fi begins to "plunge" into 
the BH with mass M. By matching to a sum of quasinor- 
mal modes (whose complex frequencies are determined by 
the mass and spin of the BH merger remnant) near the 
"mdeformed" light-ring associated with the EOB metric, 
the EOB approach produces a complete waveform that 
describes the inspiral, merger, and ringdown phases of 
binary BH coalescence. 

The EOB formalism is also highly modular: on top 
of the "base" -EOB (consisting of the 3PN EOB effective 
Hamiltonian), various elements can be added and their 
parameters adjusted. For example, a "pseudo-4PN" term 
can be added to the EOB metric, and its coefficient can 
be adjusted to help improve agreement with NR simula- 
tions. One can also add correction terms or multiplicative 
factors to the dissipative dynamics or to the amplitude 
of the waveform modes. Some of these terms attempt to 
improve agreement with the exact NR results by making 
educated guesses about some of the (non-PN-expanded) 
physics implicit in the two-body dynamics; other terms 
contain additional free parameters that can be adjusted 
to agree with NR simulations. This modularity gives the 
EOB approach a great deal of power and allows it to 
match the time-domain NR waveforms with high accu- 
racy [13, HE HO, [H-Hl] . While the EOB formalism thus 
serves as a framework to generate fast waveforms that 
agree well with NR simulations, it is less clear if the EOB 
approach — via its mapping of the two-body PN dynamics 



3 Recall that r\ = mirri2/M 2 = q/(l + q) 2 < 1/4 is the reduced 
mass ratio, and is denoted v by some authors. 



3 



onto motion in a deformed Schwarzschild background — 
provides a deeper understanding of the two-body dynam- 
ics than is already afforded by the ordinary PN equations 
of motion. This issue has been addressed in several pa- 
pers by Blanchet [l2l l6p - [63| and will be discussed in more 
detail below. 

Recently, Yunes et al. [iH (see also [HE]) have ap- 
plied the EOB formalism to model EMRIs. They show 
that if three fitting parameters are introduced into the 
dissipative portion of their model, they can accurately 
match the waveforms generated by Hughes's BH pertur- 
bation theory code [H, [67} for quasicircular inspiral in 
the Schwarzschild spacetime (with errors of < 0.1 rad 
in the phase and 0.002 in the fractional amplitude over 
a two-year integration). In retrospect, this agreement 
is not necessarily surprising because (i) in the q — > 
limit, the conservative dynamics for the EOB and BH 
perturbation theory approaches are the same — circular 
geodesies in Schwarzschild; and (ii) the dissipative dy- 
namics is solely described by the GW energy flux dE/dt, 
which is known analytically to 5.5PN order in the test- 
mass limit [68]. Combined with analytic BH absorption 
terms at 4PN order [H, [70| and adjustable parameters 
for the 6PN and 6.5PN terms in the flux, it is plausible 
to expect that such a high-order PN expansion of the 
energy flux (combined with Pade resummation and the 
factorization of an adjustable pole parameter w po io) can 
be made to match well with the numerical BH perturba- 
tion theory computation of the energy flux. Nonetheless, 
the work in Ref. [64| is an important demonstration that 
just as the EOB approach has been shown to be adept 
at fitting the results of NR simulations, it can also be 
applied to fit the results of BH perturbation theory cal- 
culations (although it remains to be seen how well the 
approach will work for eccentric, inclined Kerr orbits, 
which are expected to be typical for EMRIs). 

Yunes et al. [H, HI] also use their EMRI/EOB ap- 
proach to investigate higher-order GSF contributions 
that are not contained in Hughes's BH perturbation 
code (which only accounts for the leading-order, orbit- 
averaged dissipative piece of the GSF). In particular they 
find a phase error (between the EOB and BH perturba- 
tion theory waveforms) of ~ 3 radians due to the con- 
servative GSF terms in the EOB Hamiltonian, and < 20 
radians when the higher-order dissipative GSF terms in 
the EOB dynamics are included (see Fig. 3 of [HJ]). How- 
ever, it is far from clear that the EOB approach can make 
accurate statements about higher-order GSF effects (as 
Yunes [65[ himself notes), and part of the motivation of 
this work is to assess the degree to which the EOB for- 
malism embodies conservative GSF effects. 

By construction, the EOB formalism completely ac- 
counts for the conservative dynamics in the test-mass 
limit. The finite test-mass information it incorporates 
originates from the ordinary PN two-body dynamics, 
which is known to converge slowly in the small-mass- 
ratio limit. Therefore, even though the EOB conserva- 
tive dynamics is exact in the q — > limit, and matches 



NR calculations reasonably well in the equal-mass limit 
(with the help of extra "flexibility" parameters) , it is not 
at all clear how accurate the conservative EOB dynam- 
ics is for small (but nonzero) values of the mass ratio 
q. This issue is investigated here in the context of the 
recent GSF ISCO calculations (see also related work by 
Damour fn\). 

In particular, this study is especially interested in the 
performance of the "base" EOB 3PN Hamiltonian. As 
discussed here and in [26|, [7l[ , additional parameters can 
be introduced in the EOB formalism and calibrated to 
reproduce the results of GSF calculations. However, 
these GSF calculations are themselves currently limited 
to computing first-order in q corrections to geodesic mo- 
tion. In some situations (e.g., intermediate- mass-ratio in- 
spirals or IMRIs and possibly EMRIs) 0(q 2 ) corrections 
are expected to be important, and no "exact" numerical 
technique exists to treat this case (but see [l8l - [20j for a 
first attempt). In this case one cannot expect to be able 
to fully calibrate the EOB formalism. If one would like 
to apply the EOB approach to model IMRIs 64], it is 
important to gain insight into how the conservative EOB 
dynamics performs in the absence of any calibration to 
known numerical results. 

B. Previous self-force comparison studies 

Comparisons between PN results and BH perturbation 
theory calculations have a long history (see Refs. 16-23 
of Ref. 0]), but these involve only dissipative self- force 
effects like the radiated energy and angular momentum. 
Conservative GSF effects have been computed only very 
recently, and the first comparisons with PN results were 
performed by Detweiler [28[. In particular, Detweiler 
identified two well-defined, gauge-invariant quantified 
that could be calculated analytically in the PN approach 
and numerically in the GSF approach. These quantities 
are the angular frequency f2 of a particle on a circular 
orbit as measured by a distant observer, and the time- 
component of the particle's four- velocity u\ = dt/dri. 
(The quantity u\ can be identified with the redshift of 
a photon emitted by the particle and received by a dis- 
tant observer on the z-axis perpendicular to the circular 
orbit.) While f2 and u\ are themselves functions of gauge- 
dependent quantities (such as the orbital radius and the 
metric perturbation), one can calculate both quantities 
numerically for a particular choice of gauge in the GSF 
approach, and also analytically compute u\ as a func- 
tion of O in a PN analysis. The redshift function can be 



4 The quantities are gauge invariant in the following sense: for qua- 
sicircular orbits in Schwarzschild, the quantities and (and 
Q = u^/u\) are unchanged under an infinitesimal coordinate 
transformation x^ — > x^ + £'' provided that (i) k a d a = dt + f2c\> 
remains a helical Killing vector in the perturbed spacetime on 
a dynamical (orbital) time, and (ii) that the gauge change pre- 
serves the reflection symmetry across the equatorial plane |2ql . 
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expressed as 



-(Sa/3)l<< 



-1/2 



(1.2) 



where vf = dyf/dt = (l,v\) is the PN coordinate ve- 
locity of the particle and (g a p)i is the regularized metric 
evaluated at the particle's position. 

Using the near-zone metric previously computed to 
2PN order, Detweiler [28[ constructed the PN expansion 
of u\(y) [where y = (m2^l) 2 ^ 3 ] and compared with his 
numerical GSF evaluation. He showed good agreement 
at the 2PN level, and made a prediction for the value 
of the 3PN coefficient in u\ [Eq. 11.51 below] . In later 
work Blanchet et al. [72| extended the computation of 
(gafi)i to 3PN order and improved the accuracy of the 
GSF calculations reported in [28| , finding excellent agree- 
ment between the 3PN coefficient and a fit of its value 
to the GSF numerical results. Their refined numerical 
GSF results motivated Blanchet et al. [zl, [zH to fur- 
ther push their PN computations of (g a /3)i and u\{y) to 
higher orders. In particular they computed the logarith- 
mic corrections to the near-zone metric at 4PN and 5PN 
orders (the nonlogarithmic corrections are more difficult 
to compute and remain unknown). Their result for u\ 
(when expanded in powers of the mass ratio q) takes the 
form [Zl,l2i 

«! = 4chw + <?4f + <7 2 "PSF + 0(q 3 ), (1.3) 

where the Schwarzschild result is known exactly, 

4 chw = (1 - iyy 1 ' 2 , (1.4) 
and the leading-order self-force piece is given by 
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, (1.5) 



where the unknown coefficients a a, 015, a§, and fa at 
4PN, 5PN, and 6PN orders were determined by least- 
squares-fitting to the accurate GSF results (see Table V 
of [73[; if fa is set to zero, a value for 0:7 was also de- 
termined, but including fa caused the fits to worsen). 
Their results show that the successive PN approxima- 
tions smoothly converge to the exact GSF results (see 
Fig. 2 of Ref. [73]). The post-self- force piece Up SF was 
also calculated to 3PN order (and the logarithmic terms 
to 5PN order) , but there is no second-order GSF formal- 
ism with which to compare them. 

In addition to the comparisons of the GSF with the 
PN expansion for u\(y), Barack and Sago also computed 
u\(y) using their independent GSF code H and com- 
pared their results with Detweiler's code [3fJ. Although 



the two codes use different interpretations of the per- 
turbed motion, different gauges for the metric perturba- 
tion, and different numerical techniques, their results for 
u\(y) agree to within numerical errors. 

Recent work by Damour (7l| investigated conservative 
GSF corrections in the EOB approach. Damour points 
out that GSF calculations can provide information on 
the two functions that appear in the EOB effective met- 
ric. These functions are expanded in Taylor series in 
u = M/r (which are further resummed via Pade approx- 
imants). While pseudo-4PN and 5PN terms in this se- 
ries have been constrained by NR simulations [5(3, [Hj], 
Damour discusses how GSF calculations can similarly 
constrain higher-PN-order terms when these effective- 
metric functions are expanded in the small- 77 limit. Some 
of these constraints arise from the conservative correc- 
tions to the IS CO computed by BS (thi s is discussed fur- 
ther in Sec. IVI Al below). Damour [7l| also investigates 
how orbits with small eccentricity, as well as a special 
class of zoom-whirl orbits, can further constrain parame- 
ters appearing in small- 77 expansions of the EOB effective 
metric. The study presented here contains some overlap 
with Damour's work [Til ] . but here we focus primarily on 
comparisons of ISCO calculations with a large variety of 
PN-based methods in addition to the EOB approach. 

Even more recently, Barack, Damour, and Sago (26j 
have computed the GSF correction to the rate of pe- 
riastron advance in the small-eccentricity limit. They 
compare their numeric results with a particular gauge- 
invariant function p(x), which is related to the ratio of 
the radial and azimuthal orbital frequencies and depends 
on the small-mass-ratio expansion of functions appear- 
ing in the EOB metric. They find very good agreement 
with the known 3PN expansion of p(x), and are able to 
set constraints on higher-order nonlogarithmic terms in 
p(x) (the 4PN and 5PN logarithmic terms having been 
recently computed in [75[ and reported in [26j]). 



Previous comparisons of PN-based ISCO 
calculations with numerical relativity 



While comparisons between conservative GSF and PN 
results are very recent, comparisons between PN and NR 
results have a long history. Particularly relevant to this 
study are comparisons between PN-based ISCO calcula- 
tions and earlier work in NR involving quasicircular ini- 
tial data (QCID) calculations. By numerically construct- 
ing sequences of quasicircular initial-data sets, several 
QCID studies computed the ISCO frequency for equal- 
mass binary BHs [76H86| . These calculations typically 
involve solving a subset of the full Einstein equations 
subject to certain approximations (such as the presence 
of a helical Killing symmetry, or specifying the confor- 
mal spatial metric to be flat or a linear superposition of 
two Kerr BHs.) Several of these works also compared 
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their results with PN ISCO estimates!! The ISCO in the 
unequal-mass case has not been studied as thoroughly, 
but see Ref. [8!| for an extension of the work in [78[ to 
unequal-mass BH binaries. 

Blanchet [6l[ compared a variety of PN methods for 
computing the ISCC0 to the numerical result from [79| 
for corotating, equal-mass binaries. However, aside from 
the PN energy-function approach (to which Blanchet 61] 
derived the spin-corrections), the other ISCO estimates 
were computed only for nonspinning BHs, so it is unclear 
how to precisely evaluate all of the resulting comparisons. 
Damour et al. [§(| extended EOB computations of the 
ISCO to corotating binaries, and c omp ared their calcu- 
lations with QCID results from [jf| [la[7^| in the nonro- 
tating case. In Table HI] below I give an update of these 
equal-mass comparisons, showing how a larger variety of 
PN methods compares with more recent QCID calcula- 
tions for nonrotating BHs. My results are consistent with 
and expand on those presented in [6l|, |90| • 

In Refs. [H, |63| Blanchet reviews his results from [6l[ 
and argues against the notion that the two-body prob- 
lem in the comparable-mass limit is better represented 
by resummation methods (Pade approximants and EOB) 
than by standard Taylor PN expansions^ His argument 
is based on an estimate of the radius of convergence of 
the PN expansion of the circular-orbit energy [Eq. fl3.1[) 
below]. In the test-mass limit this radius of convergence 
is found to be x = (Mil) 2 / 3 = 1/3, corresponding to 
the frequency of the Schwarzschild light-ring (the inner- 
most location where circular orbits can exist). But in 
the equal-mass case this estimate of the convergence ra- 
dius occurs at x ~ 2.88, implying that there is no notion 
of a deformed light-ring in the comparable-mass case. 
This is in contrast to the EOB approach, which describes 
the two-body problem as containing an 77-deformed light- 
ring. Blanchet concludes that the PN two-body problem 
does not appear to be "Schwarzschild-like." Blanchet 
also argues that the 3PN value of the ISCO frequency in 
the equal-mass (rj — 1/4) limit [x; sco (l/4) 0.2, as com- 
puted from the minimum of the orbital energy] is likely to 
be accurate because the ISCO lies well within the radius 
of convergence [x(l/4) ps 2.88] of the PN series. This is 



in contrast to the test-mass limit, where the exact ISCO 
frequency xf^" = 1/6 is rather close to the light-ring 



af r chw = 1/3 



For even earlier work on the ISCO in comparable-mass binaries, 
see Refs. [87l. |88|| . Note also that Ref. [591 compared some PN 
and EOB ISCO methods, but did not include comparisons with 
numerical calculations. 

Specifically, Blanchet |6l| considered the standard PN energy 
function, EOB, the e-method, and the j'-method; these are dis- 
cussed in detail below. 

Focusing on the energy flux rather than the ISCO. Mroue et 
al. [3^1 also examined the role of Pade approximants versus Tay- 
lor expansions. They argue that Pade approximants do not al- 
ways help to accelerate the convergence of the energy flux in 
either the test-mass or equal-mass limits (although their Fig. 5 
indicates that some Pade approximants of the flux perform bet- 
ter than Taylor expansions). They also argue that the use of 
Pade approximants in generating waveform templates does not 
offer significant advantages over using Taylor expansions. 



In a more recent study of QCID, Caudill et al. [sB] im- 
prove upon the previous work of [8l| and compare their 
QCID calculations to PN formulas for the ISCO. Their 
value for the ISCO frequency for nonspinning binaries 
[A/f2i sco (l/4) s» 0.12] was found to agree more closely 
with a standard 3PN ISCO estimate j6l| [with ss 7% er- 
ror using the 3PN orbital energy; Eq. (|3.1[) below] than 
with a 3PN EOB estimate [§3] of the ISCO (with « 36% 
error; see Figs. 15-17 and Table II of 86], and Table III 
of ||). In the corotating case [MSl™ r ot (l/4) » 0.107], 
the 3PN EOB ISCO (90| performed better (ra 9% error 
versus 17% for the standard 3PN case [6l|; see Table III 
of [13]). These conclusions are qualitatively consistent 
with Fig. 3 of [13]. Comparisons between PN and QCID 
calculations of the equal-mass ISCO will be further ad- 
dressed in Sec. |W] (see Table |H] below). 



D. Summary 

In Sec. HI] I briefly review the definition of the 
Schwarzschild ISCO and the difference between the ISCO 
and the ICO (innermost circular orbit), and give a short 
discussion of the GSF. I elaborate on how the dissipa- 
tive self- force affects the ISCO, and then review the con- 
servative ISCO shift calculations by BS. I also review 
D amour's 
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reformulation of the BS ISCO frequency 
shift into the standard PN notation and gauge. 

Section HII] reviews all of the PN/EOB-based ap- 
proaches for computing the ISCO: (1) Section fill Al dis- 
cusses the minimization of the standard 3PN energy func- 
tion. Section IIIIBI examines (2) a stability analysis of 
the standard 3PN equations of motion, leading to two 
algebraic equations for the ISCO radius and frequency 
that must be solved numerically. It also discusses ad- 
ditional approaches from Blanchet and Iyer [62j. These 
involve expressing an analytic condition for the ISCO 
as a PN expansion in terms of either (3) the harmonic- 
gauge radial coordinate used in the Lagrangian form of 
the 3PN equations of motion or (4) the Arnowitt-Deser- 
Misner (ADM) radial coordinate used in the Hamiltonian 
formulation of the equations of motion. As discussed in 
[62| . these last two ISCO conditions can be reformulated 
in terms of (5) a single gauge-invariant analytic condi- 
tion that can be solved for the ISCO frequency. This 
ISCO criterion [Eq. (I3.7[) ] is particularly special because 
(i) it contains the exact Schwarzschild ISCO without ap- 
plying any resummation methods, and (ii) it produces 
the closest agreement to the BS conservative GSF ISCO 
shift of any 3PN-order method. Section fill CI computes 
the ISCO by (6) numerically solving an algebraic system 
derived from the 3PN Hamilton's equations. 

Section IIII Dl discusses a variety of "hybrid" methods 
that were originally inspired by Kidder, Will, and Wise- 
man (KWW) [9l| . These hybrid methods all involve re- 
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moving test-mass-limit terms in PN expressions and re- 
placing them with the equivalent (fully relativistic) ex- 
pressions from the Schwarzschild spacetime. In particu- 
lar, this section considers (7) a hybrid 3PN energy func- 
tion (Sec. IIII D (8) the KWW equations of motion, 
extended to 3PN order (Sec. IIII D 2"]) . and (9) a Hamilto- 
nian version of the KWW equations 92] (also extended 
to 3PN order; IIII D 31) . Section HITEl discusses two ad- 
ditional resummation approaches based on minimizing 
Pade approximants of (10) an "improved" PN energy 
function (the e-method [93[) and (11) a function based 
on the orbital angular momentum (the j -method (60j). 

Section IIII Fl discusses the EOB approach and, in 
particular, reviews the three ways of expressing the 
EOB effective-metric function A(r) (which determines 
the EOB ISCO): (12) as a Taylor series expansion, (13) 
as a Pade approximant of the Taylor series, and (14) via 
a new logarithmic expression introduced in [94j ]. Lastly, 
Sec. IIII Gl discusses (15) the Shanks transformation, a se- 
ries acceleration method that is applied to several of the 
above PN-based ISCO calculations (which are themselves 
each computed at multiple PN orders). 

Section HVl gives the results from these 15 ISCO calcu- 
lations. These are provided in Table U which shows the 
0(n) conservative correction to the exact Schwarzschild 
ISCO frequency, as well as its deviation from the exact 
BS result. Figure [T] illustrates a subset of this table in 
graphical form. The ISCO in the equal-mass case is also 
tabulated, and compared with the QCID results from 
86]. Figure [5] shows the ISCO frequency for several of 
the considered methods as a function of rj. 

Section IVl draws a number of observations from Tables 
U and HU Readers wishing to get to the main point of 
this paper as quickly as possible can skip directly to that 
section. The most important points to take away are: 

(a) The best agreement with the BS results comes 
from an EOB approach that includes a pseudo-4PN 
term that is calibrated to the comparable-mass Cal- 
tech/Cornell NR simulations in |54| . 

(b) If we do not consider calibrated methods but rely only 
on our current 3PN-level understanding of the con- 
servative two-body problem, the best agreement with 
the BS ISCO shift comes from the gauge-invariant 
ISCO condition of Blanchet and Iyer [62j]. This ISCO 
method is special because it already contains the 
exact Schwarzschild ISCO without introducing any 
"manual" resummation. 

(c) An extension of this gauge-invariant ISCO condition 
to spinning binaries shows that it also (i) reproduces 
the Kerr ISCO to the expected order in the spin pa- 
rameter, and (ii) reproduces the conservative shift in 
the ISCO due to the spin of the test-particle. This is 
discussed further in a companion paper (95|. 

(d) If we compare PN/EOB approaches with numerical 
relativity calculations of the ISCO based on the qua- 
sicircular initial-data calculations of [86], then the 



ISCO computed from the standard 3PN energy func- 
tion performs better than all EOB methods. 

(e) In both the extreme-mass-ratio and equal-mass cases, 
the 3PN EOB approach provides a single consistent 
method that computes conservative corrections to the 
ISCO in both the small-mass-ratio and equal-mass 
cases with nearly the same error (~ 27%). Calibrat- 
ing a pseudo-4PN term to NR simulations or to the 
BS result reduces this error in both limits. 

Section |VI] discusses various ways that ISCO compu- 
tations from GSF and NR can improve GW templates. 
Section IVl Al discusses how these numerical ISCO calcu- 
lations can be used to fit pseudo-4PN parameters appear- 
ing in the EOB effective metric. Section IVl Bl discusses 
how GSF results can be used to fix some of the free 
parameters that appear in phenomenological inspiral- 
merger-ringdown templates. Section IVl CI discusses how 
GSF and QCID ISCO results can constrain the undeter- 
mined functions appearing in the 4PN- and 5PN-order 
pieces of the PN orbital energy [Eq. (|3.1|)]. Section IVTDl 
briefly introduces a new approach to combine results from 
full-NR evolutions with QCID calculations to help fix 
higher-PN-order terms in inspiral templates. Section IVIII 
discusses conclusions and future work. 

Geometric units (G — c — 1) are used throughout 
this work. Note that since formulas are used from a 
large number of references, similar or identical quantities 
are sometimes denoted differently in different subsections 
of this paper. I have largely tried to keep the notation 
consistent with the original source rather than unify the 
notation throughout the paper. The context hopefully 
makes the intended meaning clear. 

II. THE ISCO AND ITS SELF-FORCE 
CORRECTIONS 

A. The Schwarzschild ISCO 

The notion of an ISCO arises from the geodesic dy- 
namics of a te st p article in the Schwarzschild geometry 
(cf. Ch. 25 of (96|). In the Schwarzschild geodesic equa- 
tions 




= E 2 -V oS (r,L), (2.1a) 



dip L dt E tn „, . 

-r = -2> i~ = t— 5 — r» (2 - lb) 

dr r z dr 1 — Zm,2/r 

the radial motion is governed by an effective potential 

V cS = (l-2m 2 /r)(l + L 2 /r 2 ), (2.2) 

where E = Ejm\ and L = L/m\ are the energy and 
angular momentum per unit test-mass, r is proper time 
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along the geodesic, and (t, r, tp) are Schwarzschild coordi- 
nates. From the condition for circular orbits, dV e g/dr — 
0, one easily finds that the angular momentum Lq and 
radius tq for circular orbits are 



m,2rl/(ro — 3m 2 ) and 



ro 



/ -■ 

2m 2 



l + Jl-12ml/L 2 



(2.3) 



(2.4) 



The latter expression indicates that there is a minimum 
angular momentum for circular orbits, L\ > L^ rit = 
12m\, and to this minimum there corresponds a radius 
ro = 6m 2 below which no stable circular orbits can exist 
(an ISCO). [Unstable circular orbits can exist below this 
radius (down to ro = 3m2), but their angular momenta 
are greater than L cr ; t .]. 

A critical radius for circular orbits can also be derived 
by finding the radius that minimizes the orbital energy of 
the test mass along a sequence of circular orbits. The en- 
ergy along circular orbits is easily found by substituting 
dr/dr = and L — Lq in Eq. (|2. la[) . yielding 



(r - 2m 2 ) 2 
r (r - 3to 2 )' 



(2.5) 



from which one can easily verify that the energy mini- 
mum occurs at ro = 6m 2 . 



B. ISCO vs ICO 

The critical point obtained by minimizing the energy 
of a circular orbit is sometimes referred to as the ICO (in- 
nermost circular orbit jf| [6l|. To clarify, an ICO is defined 
to be the frequency where the circular-orbit energy satis- 
fies dE/dVt — 0. The ISCO, on the other hand, refers to 
the point of onset of a dynamical instability in the equa- 
tions of motion for circular orbits. In Schwarzschild the 
ICO and the ISCO are clearly the same. This is also true 
for Kerr and for the conservative orbital dynamics de- 
fined via the EOB approach. In Sec. IV A 2 of Ref. (97). 
the authors show (in a Hamiltonian formulation) that the 
ISCO and ICO are formally equivalent if the Hamiltonian 
is known exactly. However, the ICO and ISCO yield dif- 
ferent critical frequencies because the Hamiltonian (or, 
equivalently, the energy or the equations of motion) are 
known only to some finite (say nth) PN order. Hence the 
analytic conditions defining the ICO and ISCO, having 
different functional forms, differ (when PN-expanded) in 



terms at n + 1 and higher PN orders. Since the condi- 
tions that define the ISCO or ICO are usually solved nu- 
merically, these implicit higher-PN terms (which would 
normally be truncated in an analytic expansion) cause 
the numeric values for the ISCO and ICO frequencies to 
differ. (See also Sec. 4.3 of [§8[ for a related discussion.) 

From a practical standpoint, both an ICO and ISCO 
signify a frequency at which the standard PN descrip- 
tion of adiabatic circular orbits breaks down. In studies 
that consider the finite- mass-ratio case (see e.g., the ref- 
erences in Sec. II CI above') . several papers (e.g., foil l90l|) 
take the viewpoint that the "correct" quantity to com- 
pare with comparable-mass QCID simulations is the ICO 
rather than the ISCO (although the QCID papers 
|86| themselves usually refer to this quantity as an ISCO). 
This is in part because an ISCO (in some approaches) is 
not always defined at each PN order or for 77 ~ 1/4 (see, 
e.g., Sec. Hill below or Sec. IV A 2 of [97]). However, in 
the small-ry limit there are both ICO and ISCO methods 
that yield well-defined (but not always well-behaved) re- 
sults. This will be made clear in Sec. IIVI below. In the 
rest of this paper, for simplicity I will refer to both ICOs 
and ISCOs as "ISCOs." The context will make clear if 
the method in question is formally an ICO or an ISCO. 



A (very) brief overview of the gravitational 
self-force 



How does the ISCO change when the mass mi is no 
longer completely negligible? The answer to this question 
is the purview of self-force calculations. The self-force 
causes the motion of the point-mass mi to deviate from 
geodesic motion via 



dT 2 



dx@ dx 1 
dr dr 



clf(l) 



- t sclf(2) 



-0(q 3 ), (2.6) 



where on the right-hand-side the self-force is expanded in 
powers of the mass ratio [a" cl ^ n ^ cx q n ]. The background 

metric used in the left-hand-side is usually taken to 
be the Schwarzschild or Kerr metric. The full spacetime 
metric includes perturbative corrections of the form 



9V 



where hffl oc q n . 



0(q 3 



(2.7) 



The leadin g-or der GS F ha s been derived by several 
authors (see |2ll424l l99l . 1 1 OOj ] for references and recent 
work) and is given by 



<lf(i)(*) = V a ^ 7 (*), 



(2.8) 



In Ref. [93 the ICO is referred to as the MECO (maximum 
binding-energy circular orbit). The designation LSO (last stable 
orbit) is used by some authors to refer to an ICO, and by others 
to refer to the generalization of the ISCO to generic orbits. 



where V a ^ 7 is a differential operator proportional to a 
covariant derivative, and (z) is the regularized met- 
ric perturbation evaluated at the position z M of mi (an 
overbar means to take the trace-reversed part). This reg- 
ularized metric perturbation is a smooth solution of the 
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homogenous linearized perturbation equations. It is con- 
structed from the first-order retarded metric perturba- 
tion h™* by analytically subtracting out a singular con- 
tribution /i| . The retarded metric perturbation is itself 
a numeric solution of the inhomogeneous linearized per- 
turbation equations with a point-particle source. Note 
that the motion of mi is equivalently described by purely 
geodesic motion, but in terms of a background metric 
g° + and a new proper time r', 

d 2 x a „ . n p. dx^ dx" 1 , N 

l^+^° + ^] — —=0. (2.9) 

For further details see [IH [24[ and references therein. 

Self-force effects are more easily studied by splitting 
the GSF into a dissipative (time-odd) and conservative 
(time-even) piece. If one views the GSF as moving the 
particle m\ along a sequence of geodesies instantaneously 
tangent to its motion, then the dissipative and conser- 
vative pieces of the GSF modify the constants of motion 
parametrizing these geodesies. The dissipative GSF piece 
is responsible for secular changes in the "intrinsic" con- 
stants of the motion: the energy E, azimuthal angular 
momentum L, and the Carter constant Q. These changes 
give rise to radiation-reaction, causing the orbital sepa- 
ration, eccentricity, and inclination to slowly evolve on 
a radiation-reaction time scale. The conservative GSF 
also modifies the intrinsic constants of the motion, but 
does so in an oscillatory manner that averages to zero 
on an orbital time scale. However, the conservative GSF 
can also affect the "extrinsic" constants of the motion: 
these constants are responsible for the orientation of the 
geodesic and the location of the particle on the geodesic. 



D. Dissipative self-force effects on the ISCO 

The effect of the dissipative GSF on the ISCO was ad- 
dressed in a study by Ori and Thorne [lOl| ] . Specifically, 
they showed that the region near the ISCO gets "blurred" 
into a transition regime lying between the adiabatic in- 
spiral and plunging phases. They derive approximate 
equations of motion for the transition regime by expand- 
ing the geodesic equations about small deviations from 
the ISCO. Dissipative GSF effects are included by allow- 
ing for dissipative changes in the E and L that enter the 
effective potential V e s (note that for Kerr, V e s depends 
on both L and E). They find that the radius of the 
transition regime Ar = rj sco — rt rans and the shift in the 
or bital frequency Afl = fitrans — Slisco is given by (Sec. IV 
of flOlj) 

A " S 4.4g 2 / 5 . (2.10) 



with the corresponding fluxes at infinity given by 



— 18g 2/5 and 
m 2 



n,, 



Note also that the energy and angular momentum 
radiated during the transition regime is given by 
[cf. Eq. (3.26) of pl| ] 



— « 0.096 9 9 / 5 

TO2 



and 



ml 



(2.11) 



AE , 
— — w 0.000 97q 2 
At H 



and 



AL/m 2 
At 



O.OUq 2 . (2.12) 



In Eqs. (|2.10|) - (|2.12j) . the numerical prefactors assumed 
the Schwarzschild spacetime; the corresponding values 
for equatorial orbits in Kerr are easily derived from [ 1 ll j . 
These scalings for the transition region were also inde- 
pendentl y d erived in the EOB analysis of Buonanno and 
Damour [57| (which considered the nonspinning case but 
for arbitrary mass ratios). 



E. The Barack-Sago conservative GSF ISCO shift 

To compute the conservative GSF corrections to the 
ISCO, Barack and Sago [H, 0] analyzed the equations of 
motion in the form 



d r 
df 2 
dE 
~dJ 



ldV oS (f,L) 



Of 



dL 



(2.13a) 
(2.13b) 



where the hats refer to quantities along the new orbit of 
the particle (which is no longer a geodesic and on which 

the energy and angular momentum parameters E and L 
are no longer constants of the motion). These equations 
are then expanded in terms of a slightly eccentric orbit 
near the Schwarzschild ISCO. The resulting shifts in the 
ISCO radius and orbital frequency are then expressed 
in terms of the components of the GSF evaluated at the 
ISCO. The difficult part of the BS analysis is numerically 
computing the components of the GSF along circular and 
eccentric geodesies (see Eh3 for the details) . Working in 
Lorenz gauge, BS fincfl [l|, bl 



6m 2 - 3.269(±0.002)mi, 



(2.14a) 



™ QLorcnz 
m 2"isco 



6- 3 / 2 [l + 0.4869(±0.0004)g]. (2.14b) 



As pointed out by Damour [7l| (see also Sec. Ill D of 
[l02| and Sec. Ill B of [30]), the Lorenz gauge used for 
the calculations in BS is not asymptotically flat. Rather 
the asymptotically-flat time coordinate is related to the 
r — > oo limit of the perturbed binary metric by cftfiat = 
(1 + 2a) 1/2 dt Lorcnz , where 



m\E\ 



r (l - 2m 2 /r ) 



(2.15) 



9 Recall that most of the self-force literature uses (/i, M) in place 
of (mi , m.2). Because we compare with PN computations, we use 
here the conventions commonly employed in the PN literature: 
M = mi + r«2, fi = mxmi/M, r\ = fi/M, and q = mi/m,2 < 1. 
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Here ro refers to the Schwarzschild coordinate radius of 
the particle's (mi) circular orbit, and 



Ei 



Ei 



1 - 2m 2 /r 



mi y/l - 3m 2 /r 



(2.16) 



is the particle's conserved energy per unit rest mass. 
Evaluated at the ISCO (ro — 67712), a takes the value 
q^/2/6. To convert angular circular-orbit frequencies 
from Lorenz coordinates (S^Lorcnz = dip / dt-L, oleaz ) to 
asymptotically- flat coordinates (Ofl a t = dtp / dtfat) , we 
use [nl 



Qflat = ^ Lmc " z fa ft Lorcnz [l - a + 0(q 2 ) 
V 1 + 2a 



(2.17) 



When converted to asymptotically-flat coordinates, the 
frequency of the ISCO becomes 



Oisco.sf 
- flat 



g-3/2 



1 



„BS 



(2.18) 

where c^ s = 0.4869 (±0.0004) Q. To simplify compar- 
isons with PN ISCO calculations, it is convenient to ex- 
press this result in terms of M and 77 instead of 777.2 and 
q. Multiplying both sides of Eq. (j2"TT8"| by M /m 2 and 
using 77 = q + 0(q 2 ), we have 



Mil 



where 



flat 



1 



6-3/2 [l + c ^ + 0(77% (2.19) 



„BS 



1.2512(±0.0004). (2.20) 



It is this number that will be compared with the PN- 
based ISCO calculations below. 

Note that the ratio of the conservative to the dissi- 
pative [Eq. (I2.10[) ] changes to the ISCO radius or fre- 
quency is roughly equal to ~ O.lg 3 / 5 . For small mass 
ratios (q ~ 10~ 2 -10~ 7 ) this implies that the "blurring" 
of the ISCO by the dissipative GSF overwhelms the 
small conservative-GSF shift in the ISCO by a factor 
of - 10 2 -10 5 . While the dissipative GSF ISCO shift 
becomes more important than the conservative shift as 
the mass ratio gets smaller, in the comparable-mass limit 
the entire notion of an ISCO is not well-defined (at least 
in the presence of dissipation). The conservative-GSF 
ISCO-shift is therefore unlikely to be a quantity of ob- 
servational importance. Rather its importance lies in the 
fact that it represents a unique, strong-field critical point 
in the conservative two-body dynamics that can serve as 
a test of numeric GSF codes and a point of comparison 
with PN (and perhaps NR) calculations. 



III. A CATALOG OF PN-BASED METHODS 
FOR COMPUTING THE ISCO 



a Maple code was developed to numerically compute 
the ISCO frequency £lisco(v) as a function of the reduced 
mass ratio. 



A. PN energy function 

One of the simplest methods for determining the ISCO 
(in this case an ICO) is to minimize the PN circular-orbit 
energy E PN (il) with respect to frequency. For nonspin- 
ning binaries, this energy is given by [see Eq. (3) of [6l[ 
and references therein^ 



E FN {Q) _ 1 

T)M 



x < 1 + x 



27 19 



3 T) 

"4 ~ 12 
2 
T 

24 



675 
~64~ 



34 445 205 
~96~ 
3969 



576 



7T 2 I 77 



155 
"96 



n 



35 

5184 



17 



45 927 
512 



128 

»?e 5 (»?) 



f v 448 , 
f?e 4 (7j) + —77 In x 



4988 1904 



35 



15 



-n 



77 In a; 



(3.1) 

where x = (Mfl) 2 ^ 3 and the known v alue for t he regu- 
larization parameter A = —1987/3080 [l 1 lMl 1 5j j is used 
throughout this article. Equation (|3.1j) also includes 
newly computed logarithmic correction terms at 4PN and 
5PN orders [73|. The test-mass pieces of these 4PN and 
5PN terms are known from the exact Schwarzschild ex- 
pression [Eq. (|3.2j) below]. The functions 64(77) and e 5 (?7) 
denote some unknown polynomials in 77. Section [VIC] dis- 
cusses how we can use our present knowledge of the ISCO 
from GSF and QCID calculations to help constrain these 
functions. However, the ISCO comparisons discussed in 
Sees. IIVI and fVl below will only make use of the circular- 
orbit energy to 3PN order. 

Note that in the small-mass-ratio limit, it is well- 
known that standard Taylor PN expansions converge 
slowly 043 • 

For example, Taylor expanding the 
Schwarzschild circular-orbit energy 



E Schw {n) (l-2x) 



rjM 



(1 - 3x)!/2 



1, 



(3.2) 



and computing the ISCO frequency at each PN order, one 
needs to go to at least 4PN order to reproduce the exact 
test-mass ISCO (Affi S chw = 6~ 3/2 « 0.0680) to within 
12% (see also Sec. II of Ref. [6(|). Since the (nonre- 
summed) 3PN energy function poorly matches the test- 
mass ISCO frequency [Eq. flSU predicts Mfif™^ = 
0) « 0.0867], this method cannot be straightforwardly 
compared with the BS conservative GSF ISCO shift. 



This section reviews nearly all PN-based methods for 
computing the ISCO. For each method discussed below, 



For extensions to the case of spinning binaries, see [6ll.l9Cl.H03l - 
1 106)1 . For eccentric or tidally distorted binaries, see Tl07i 4l 10(1 . 



10 



B. Stability analysis of the PN equations of motion 

The ISCO can be determined by directly analyzing the 
conservative-pieces of the PN two-body equations of mo- 
tion. In harmonic coordinates the relative two-body ac- 
celeration can be written in the form 



a = -(M/r 2 )[A(r, r, ip)n + B(r, r, <p)v] 



where M = mi + 771.2 is the total mass, r is the relative 
orbital separation, n is a unit vector that points along 
the relative separation vector, and v is the relative orbital 
velocity. The orbital phase angle is denoted tp and an 
overdot refers to a derivative with respect to coordinate 
time t. The functions A = 1 + A 1PN + A 2PN + A 25PN + 

j4 3PN + A 3.5PN and g = glPN + g2PN + g2.5PN + g3PN + 

g3.5PN are k nown to 3.5PN order (see @ for references). 
For illustration, the 1PN pieces are 

A 1PN = -2(2 + r,) — + (1 + 3r,)v 2 - ^f 2 , (3.4a) 
r I 



B 



1PN 



-2(2 - V )r, 



(3.4b) 



where v 



+ r ip for planar motion. The remaining 



pieces can be found in Eqs. (181)— (186) of @. Note that 
at 3PN order I use the form given in Eqs. (185)-(186) 
in which a gauge transformation has been applied to re- 
move the logarithmic terms. Since we are concerned only 
with conservative effects, the radiation-reaction pieces at 
2.5PN and 3.5PN orders are set to zero. 

To compute the ISCO we follow the prescription given 
in Sec. Ill A of Ref. [§lj|. This involves (i) expressing 
Eq. p. 31) as a system of first-order equations for the vari- 
ables r, lo = (p, u = f, (ii) linearizing those equations 
about a circular-orbit solution (r = u = Cj = 0), and (iii) 
finding a criteria for the stability of the solution. This 
produces a system of equations for the ISCO radius and 
orbital frequency (r ,u;o) [Eqs. (3.3) and (3.6) of [9l|] : 



w 2 



MA /rl 



w 



dA 



Aq \8uj 
M fdB 



Cn = l-2=M^ +^ ^ 



+ 2 



du 



ro_ (dA 
A \dr /tl 
1 ujq (dA 
~ 2 ~A~ Vduj 



(3.5a) 



0, (3.5b) 



where the subscript means to evaluate the quantities 
along a circular orbit. 

Solving these equations numerically yields well-defined 
solutions at 2PN order, but no physical solutions at 1PN 
or 3PN orders. Even at 2PN order, in the test mass limit 
the ISCO is found to be at r = 6.505M (see Table II of 
[9lj ) instead of the correct harmonic coordinate radius 
of 5M. Strangely, as the reduced mass ratio is increased 
from to 1/4, both the ISCO radius and the ISCO an- 
gular frequency increase. This approach does not appear 
to yield sensible results for the ISCO. 



Blanchet and Iyer [6 
approach (see also [116 [ 



pursued a variation of the above 
). Instead of solving Eqs. (13. 5[) 



numerically, they derived a PN series expansion for ujq 
and Co [a quantity equivalent to Eq. (|3.5bl) ] in terms of 
the harmonic radial coordinate^, 



2 M 

W 0,har. — ~3 1 ^ 



(3.3) + -3- 



M 3 

T 

o 



'o 
-10 



. M 2 ( 41 
3 + 7 l) + —jT 6 + —7] 



-75 707 41 2 „„, fr 

1 vr 2 + 22 In -f 

840 64 \ r' n 



+ ™2 +ri *) + 0(8)}, (3.6a) 



^har. _ M 
' n 



— ( 

ro 



M 3 



r, 



(i 



-70 



9 4 

29 927 
840 



s M 2 



451 
"64" 



30 + ^77 + rf 



22 In I ^ 



iyr, 2 -^; =i). C-S.f.l M 



Reference [62j also derived relationships equivalent to 
Eqs. (|3.6p but starting from the 3PN Hamiltonian in 
ADM coordinates [see their Eqs. (6.24) and (6.38)]. 
These Hamiltonian-based expressions uj 2 adm and Co I}M 
are expressed in terms of the ADM radial coordinate R 
and differ from Eqs. (|3 .6[) at 2PN and higher orders. 
However, Ref. [62J showed that the two sets of equa- 
tions agree if one applies the coordinate transformation 
between the ADM and harmonic coordinate radii. In 
either formulation one can solve the equation for Co to 
determine the ISCO radius and substitute the result into 
the equation for w 2 to find the corresponding orbital fre- 
quency. However, because of the coordinate-dependent 
nature of these two formulations, they yield different nu- 
merical results at 3PN ordeim: in the test-mass limit 



„3PN 



5.93M (Mwg^ r « 0.0544) while in the Hamil- 
tonian formulation i?jj PN « 5.76M (Mwf A N DM 0.0561). 
This ~ 3% difference in the frequencies presumably arises 
from differences at 4PN and higher orders. 



1. Gauge invariant description of the PN ISCO 

Blanchet and Iyer 62] also derived a gauge-invariant 
form for the ISCO condition (% = 0. This comes from 
solving Eq. (|3.6a[) for ro (or Ro in the ADM case) in 
terms of the PN frequency parameter xo = (Mwq) 2 / 3 , 



11 Note that in their derivation of Eqs. GSM , Ref. [H] used the 3PN 
equations of motion in a gauge in which the logarithmic terms 
are still present. These logarithmic terms depend on an arbitrary 
gauge constant r' associated with the choice of coordinates. Also 
note that we introduce here the notation 0(n) to refer to terms 
of order n/2 PN [i.e., 0(c" n )]. 

12 At 1PN order both methods are identical; at 2PN order there is 
no ISCO in either formulation. 
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and substituting the resulting expression for r (x) [or 
Ro(x)\ into the expression for Cft ar - (or C^ DM ). The 
result is an expression for C = C M 2 /xq [their Eq. (6.1)] 
that depends directly on the ISCO orbital frequency and 
not on any gauge- dependent radius, 



Co — 1 - 6xo 

V397 123 . 
\2 16 W 



i] — 14ri z 



x 3 + O(8). (3.7) 



This expression has the very interesting property that it 
yields the exact test-mass Schwarzschild result (cc ISCO = 
1/6) at all PN orders without any form of "resumma- 
tion. " Finite-mass-ratio effects do not enter until 2PN 
order. In this case the ISCO exists for all r\ and can be 
solved analytically to yield 



Mlo 



2PN 
ISCO 




777/12- 



14 1 
-0{v 2 )]- 



-i 3/2 



(3.8) 
(3.9) 



At 3PN order the ISCO only exists for r) < 0.183; for 
larger mass ratios all circular orbits are stable. Equation 
(I3.7[) can be solved exactly at 3PN order, but not in a 
simple form. Its expansion for small r\ is given by 



r 3PN 
l ISCO 



e-^i + c^ + otf)], 



M, > 3PN 



-3/2 



[1 



c c o3 P Nr? . 



0(ry 2 



where 



„C 3PN 



„C 3PN 



565 
432 
565 
288 



41tt 2 
1152 

417T 2 

768 



0.956 608 407. 
1.434 912 612. 



(3.10a) 
(3.10b) 

(3.11) 
(3.12) 



Note, in particular, that the coefficient c^° 3PN ~ 1.435 
differs from the exact value cff" = 1.251 by 14.7%. As 
we will see in Sees. IIVI and [V] below, this agreement is 
better than any 3PN order estimate, including all hybrid, 
resummed, or EOB methods. The extension of the above 
ISCO condition to spinning systems is discussed in [95j 
and briefly in Sec. fVl below. 



C. PN Hamiltonian 

While Ref. (HJ determined an analytic condition for 
the ISCO using the PN Hamiltonian in ADM coordi- 
nates, an alternative method follows from the work of 
Ref. [92[. Here one starts with the reduced PN Hamilto- 
nian, 



u 



u 



ADM 



Hn+Hipn+'H 



(-2PN 



u 



3PN, 



(3.13) 



where the Newtonian and 1PN terms are 

p2 X 



n 



1PN 



1 1 

2~I 



[(3 



v)P 2 



R 



and 



(3.14a) 



1 1 



(3.14b) 



R = RM is the ADM radial coordinate, P is the con- 
jugate momenta divided by the reduced mass /i, and 
TV = R/ R is the unit orbital separation vector. The 
2PN and 3PN terms can be read from Eq. (5.9) of [62 1 
or references therein. Hamilton's equations are 



(3.15) 
(3.16) 



dR 


dH ABM 




dn ADM 


It ~ 


dP R ' 


~dt 




dP R _ 


dH ABM 


dP^, 


= 0, 


dt 


dR ' 


dt 



where (R, \1/) are ADM polar coordinates. 
The conditions for the ISCO are 



0, Pfi 



on 



ADM 



0. 



dPn 



d 2 n 



ADM 



0. 



dR "'OR dR 2 

(3.17) 

Following [92[ , the conserved energy E = [iE and angular 
momentum J = Pq, = J/iM are 



E = U{R, P) 



J = Rx P, 



(3.18) 



and the above conditions for the ISCO are equivalent to 
the system 



dE 
dR 



(R, J) = 



d 2 E 
dR 2 



(R, J) = 0, 



(3.19) 



where Eq(R, J) is the energy evaluated along a circular 
orbit and is obtained by substituting N ■ P = and 
P = J 2 /R 2 (with J = |J|) into Eq. ([3~T3]) [see Eq. (5) of 
[92]|]. For example, up to 1PN order the energy is [93] 



E (R, J) 




V 



I (sC-aC) +0(4) 
8 I i? 4 R? 



+ 0(4) 
(3.20) 



where the //-dependent terms have been separated to il- 
lustrate the hybrid method discussed below. 

By solving Eqs. (|3.19[) numerically the ISCO values 
Rq and Jo are determined. The corresponding ISCO fre- 
quency is found from 



MQ = M 



d^> 
~dt 



dE 
~dl 



(3.21) 
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At 2PN order no solutions for the ISCO were found. At 
1PN order there are solutions, although in the test-mass 
limit they differ significantly (MfiJ PN w 0.0234, 7?g PN = 
11.1) from the exact result (as expected). At 3PN or- 
der the results are somewhat better (Mf2g PN s» 0.0396, 
Rf N = 7.55) but still differ by 42% from the exact test- 
mass frequency. This method does not appear to be well- 
suited to finding the ISCO. 



D. Hybrid methods 

With the exception of the gauge-invariant method dis- 
cussed in Sec. IIIIB 11 the above methods cannot repro- 
duce the correct ISCO frequency in the test mass limit. 
This is not necessarily surprising since (i) the ISCO oc- 
curs in the strong field where the PN expansion starts to 
break down, and (ii) the PN expansion is known to con- 
verge more slowly in the test-mass limit. To help remedy 
this problem, a hybrid a pproach was introduced by Kid- 
der, Will, and Wiseman [9ll Ill7j | to enforce the standard 
Schwarzschild dynamics in the test-mass limit. The ba- 
sic philosophy behind the hybrid approach is to replace 
the test-mass limit parts of some PN expanded function 
(i.e., the leading-order terms in an expansion in rf) with 
the equivalent terms from the exact Schwarzschild repre- 
sentation of the same function. Many of the previously 
discussed methods for computing the ISCO have hybrid 
analogs which we now describe. 



1. Hybrid energy junction 

A hybrid version of the 3PN energy function in 
Eq. p.ip is easily computed by removing the test-mass 
pieces and replacing them with Eq. (|3.2[) . [Note that the 
PN expansion of Eq. Q3.2p coincides with the test-mass 
limit of Eq. (J37TJ) .] The result is 



pPN 
hybrid 

T)M 



(1 - 2x) 



(1- 

X 

~ 2 
34 445 
576 



Zx) 1 / 2 
V 

12* 
205^ 
"96"' 



- 1 



19 



155 
"96"' 



24 



35 
5184' 



(3.22) 



This is equivalent to expressions found in Ref s. 11181 . 1 1 1 Qj | 
(which were inspired by the approach in [9lT Ill7j |). Un- 
like some of the other methods investigated here, this 
hybrid energy function produces an ISCO (more appro- 
priately an ICO) that is uniquely defined at all PN orders 
and for all values of rj. Like all of the remaining hybrid, 
resummation, and EOB methods discussed below, it re- 
produces the exact Schwarzschild result for the ISCO in 
the test-mass limit. 



2. Hybrid PN equations of motion 

In the original KWW hybrid approach [9J li~i~7j , the 
conservative parts of A and B in Eq. (|3.3I) were split into 
test-mass and non-test-mass pieces. The PN test-mass 
pieces were then replaced by the exact test-mass terms 
that arise from the geodesic equation of the Schwarzschild 
metric written in harmonic coordinates (see Sec. II A of 
91]). This defines new (hybrid) equations of motion 

a H = -(M/r 2 )[{A s + A n )n + (B s + B n )v], (3.23) 

where the Schwarzschild pieces are [Eqs. (2.6) of f9lj ]] 



-4 S 



B s = 



1 - M/r 
(1 + M/r) 3 

1 - 2M/r 



M/r 



1 - [M/rf 



M 



1 - (M/r) 2 



(3.24a) 
(3.24b) 



and the non-test-mass PN pieces A v and B v are found 
by removing the ^-independent terms from A 1PN , A 2PN , 
A 3PN , B 1PN , B 2PN , and B 3PN . (As in Sec. 11111 we ig- 
nore the dissipative terms at 2.5PN and 3.5PN orders.) 
The ISCO radius and angular orbital frequency is then 
computed by numerically solving Eqs. (I3.5[) (substitut- 
ing A -> A s + Arj and B -> B s + B v ). Like all of the 
hybrid and other remaining methods below, the KWW 
approach yields the exact Schwarzschild ISCO in the test- 
mass limit. 

This hybrid approach was criticized by Ref . [92j . They 
developed a Hamiltonian formulation of the hybrid ap- 
proach (discussed in Sec. IIIID 31 below) which gives dif- 



ferent results for the ISCOJB More troubling, Ref. [91| 
found that while at 1PN order the ISCO radius moves 
inward as ry is increased, at 2PN order it moves outward 
(see Fig. 4 of Ref. HU), in contradiction to the exact 
result of Ref. [l[ . Here I extended the KWW hybrid ap- 
proach to include 3PN order terms. At 1PN and 3PN 
orders, the ISCO radius moves inward (as expected) for 
small r]; but unlike at 2PN order, as r\ increases the so- 
lutions exhibit a discontinuous jump, after which they 
move to larger radii. Above some critical value for 77, no 
solutions for the ISCO are found at 1PN and 3PN orders. 



13 The KWW hybrid approach was also criticized in Ref. |93j , which 
argued that some of the finite-?? terms in the 2PN equations of 
motion [i.e., terms in A v and B v ] amount to very large fractional 
corrections to the test-mass terms. However, this argument is not 
entirely justified. It is more appropriate to compare the ratios 
\A V /A S I or \B V /B S | , which I find are typically < 0.3 at the ISCO. 
The Schwarzschild terms do in fact dominate the finite-77 terms, 
although not by a large factor. In comparison, the 0(rr) terms 
in the EOB effective metric potential A(r) [Eq. 3-3 1 ft below] 
are much smaller than the leading-order Schwarzschild term (by 
factors < 0.009 near the ISCO) and decrease very quickly with 
increasing radius. This suggests that the EOB approach is a 
much better perturbative scheme than the KWW equations. 
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3. Hybrid PN Hamiltonian 

Another hybrid approach (suggested in [92]) involves 
replacing the test-mass pieces of the PN Hamiltonian 
with the Hamiltonian of the Schwarzschild spacetime in 
isotropic coordinates. As implemented in [921 ] . this in- 
volves replacing the test-mass terms in the first line of 
Eq. (|3.20j) with the reduced circular-orbit energy of a 
particle in the Schwarzschild spacetime [Eq. (8) of [U 



7TlSchw / T> T\ 

A),iso {H,J) — 



1 - l/(2fl) 
1 + l/(2i?) 



1+ M+ * " 



-1. 



. , .. 1/2 

2RJ R 2 ^ 

J (3.25) 
substituted into 



The resulting energy function is 
Eqs. (|3~T9l> and (j3~2Tj) to determine the ISCO. 

As pointed out in [92f . this approach yields different 
numerical results from [9l| (even after correcting for the 
change in coordinate systems). However, it does be- 
have qualitatively similar to the KWW hybrid approach: 
the ISCO radius moves inward at 1PN and 3PN orders, 
and outward at 2PN order (and solutions were again not 
found above some critical r\ at 1PN and 3PN orders). 



E. Resummation approaches: the e- and j-methods 

Additional methods for computing the ISCO — based 
on minimizing the Pade approximant J^l of certain 
functions — were introduced in Refs. [60, |93|. The e- 
method [§3| involves minimizing the Pade approximant 
of a new energy function e{x) defined by 



e(x) 



(M + E) 2 -m 2 - ml 



2mim2 



1, 



(3.26) 



where E is the orbital energy whose PN expansion for 
circular binaries is given by Eq. (I3.ip . The justifica- 
tion for this expression is discussed in [6(1 193| • The ex- 
plicit form of e(x) that is minimized to determine the 
ISCO is obtained by: (i) substituting the PN expansion 
E = E PN from Eq. ([371]) into Eq. (|3~26)) . (ii) Taylor ex- 
panding the resulting expression for e{x) in terms of x 
to 2PN order (denoting the result as T2[e(x)]) or to 3PN 
order (T^[e(x)\), and (iii) computing the Pade approxi- 
mant of the result (denoted as ep 2 (x) — Pi p2[e(x)]] or 
ep 3 (x) = P 2 1 p3[ e ( a:; )]])- Explicit expressions for ep 2 and 
ep 3 are given in Eqs. (4.6) of j6(j. They reproduce the 
Schwarzschild ISCO in the test-mass limit. (See [37j for 
a criticism of this approach.) 

In the j -method, Ref. (6(i| proposed another function — 
j(x) = J/(fiM), where J is the magnitude of the orbital 



angular momentum of the system — whose extremum also 
defines an ISCO. This function is computed by taking the 
1PN, 2PN, and 3PN Taylor series for j 2 (x) (Ti[j 2 (a;)], 
^2[j 2 (x)], T^\j 2 (x)]) and constructing the Pade approx- 
imants 3 \(x) = PfPlL^W 3%{ x ) = ^iPMj 2 (*)]], 
and jp 3 (x) = Po\Tz\j 2 (x)]]. Explicit exp ressions are 
found in Eqs. (4.16) of (gfj. Reference [6TJ argued that 
the j-method is preferable to the e-method because un- 
like the 1PN Pade approximant e Pl (x) = P 1 °[Ti[e(a;)]], 
the test-mass limit of the 1PN Pade approximant j Pi (x) 
already reproduces the Schwarzschild ISCO. Additional 
desirable properties of the j-mcthod are discussed in 
601 B 



F. EOB methods 

The effective-one-body (EOB) approach models the 
conservative two-body dynamics in terms of the dynam- 
ics of a single particle in the background of a deformed 
Schwarzschild geometry. The dissipative dynamics is in- 
corporated by supplementing Hamilton's equations with 
radiation-reaction forces [these are based on vario us wa ys 
of "resumming" the energy flux (see e.g., [541 [93l [l20j )]. 
Since we are concerned with purely conservative correc- 
tions to the ISCO, we need only consider the conserva- 
tive EOB dynamics, which satisfy Hamilton's equations 
[Eqs. (2.7)-(2.10) of [13], 



dR _ dH TeaX 
dT ~ dP R 
dip <9£P eal 



dT dP L 

dP R _ 
dT 



dH real 



dR 



(R,P Rl P v ), 
(R,P R ,P V ), 

(R, PRl Pip): 



dPp 
dT 



= 0. 



(3.27a) 
(3.27b) 

(3.27c) 
(3.27d) 



Here motion is restricted to the equatorial plane, and 
i/ rcal = [iH real is the nonspinning real EOB Hamiltonian 
[Eq. (2.11) of 13], 



H^(R,P R ,P V ) = Af^/l + 2 ?/ (-ff cff /M-l)- (3-28) 

The effective EOB Hamiltonian H cS = nH cS is [e.g., 
Eq. (5) of [m|] 



H 



off 



A(r) 1 + ^ + 



B(r) 



Pi 



(3.29) 



The Pade approximant of some function / whose power series is 
f( z ) = "^2 a nZ n consists of a ratio of power series Pp[f(x)] = 
(Ei ^» zl )/(5Zj c j zJ )■ Pade approximants are useful because 
they tend to accelerate the convergence of a series. 



Note that [60ll also defines a third invariant function for comput- 
ing the ISCO (the k-method) that is related to the periastron 
advance rate. Reference 16011 considers this method less prefer- 
able than the others so I do not consider it here. 
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with Z3 = 277(4 - 377), p v = P v /(fj,M), p r = P r /fi, 
r = R/M, and t = T/M. The functions A(r) and B(r) 
appear in the EOB effect ive metric [in Schwarzschild 
gauge; see, e.g., Eq. (1) of [HJ], 

ds 2 cS = -A(r)dP + B(r)dr 2 +r 2 (de 2 + sin 2 6dip 2 ). (3.30) 

The Taylor expansions of the functions appearing in 
the EOB metric have been computed to 3PN order and 
are given by [56|, [6(| 



and 



B{r)A(r) = D{r) 
where 

04(17) 



^ + 2(3,-26)4, 



94 41 



32 



(3.31) 



(3.32) 



(3.33) 



Note that the 1PN contribution to A(r) is exactly zero. 
We have also included a pseudo-4PN contribution to 
A{r), where the coefficient is parametrized as |54( 



05(17) = f?(Ao + A1J7). 



(3.34) 



For our calculation of the EOB ISCO we shall not need 
to make use of the functions B(r) or D(r). 

We further define the functions A T2PN (r), A T3PN (r), 
and A T4PN (r) as the Taylor expansion (13.31ft truncated 
at 2PN, 3PN, or 4PN order. We also define the follow- 
ing Pade approximants 60] to A{r) which are listed in 
Eqs. (50)-(56) of Ref. HfT 



A 



P2PN 



Pl[A 



2PNi 



r(-4 + 2r + rj) 
2r 2 + 2i] + rrj 



at 2PN order, 
(3.35a) 



A P3PN = P 3 1 [A 3PN ] = 



Num(^) 
Den(Ai) 



at 3PN order, with 
(3.35b) 



Num(^J) = r 2 [a A {ri) + 877 - 16 + r(8 - 277)], (3.35c) 

Den(^J) = r 3 (8 - 2ry) + r 2 [04(77) + 4ry] 

+ r[2a 4 (?7) + 877] + 4[t7 2 + 04(77)], and (3.35d) 



A P4PN ee P^A 4 ™] = Num( f 4) at 4PN order, with 
Den(^) 

(3.35e) 



Num(4) = r 3 [32 - 24?7 - 4o 4 (?7) - 05(77)] 

+ r 4 [04(77) - 16 + 877], and (3.35f) 



Den(^) = -04(77) - 805(17) - 804(77)77 + 205(77)77 

— I677 2 + r[— 804(77) — 405(77) — 204(77)77 — I677 2 ] 

+ r 2 [-4o 4 (7/) - 205(77) - I677] 

+ r 3 [-204(77) - 05(77) - 877] + r 4 [-16 + 04(77) + 877]. 

(3.35g) 

Note that the Taylor expansion in u — M/r of the 
above Pade approximants reduces to Eq. (13.31)) at the 
appropriate PN order. However, the Pade approximants 
of Air) also have the following interesting property: if 
one takes any of Eqs. (|3.35|) and computes the Taylor 
expansion not in it but in 77, one still arrives at the Taylor 
expansion of A(r) in Eq. (|3.31[) . 

The EOB ISCO is an inflection point in the radial mo- 
tion given by [Eq. (2.17) of (13]] 



dH 



:cal 



2 trrcal 



-(R,P R = 0,J) = = 



d 2 H 



■(R,P R = 0,J), 



dR v ' x ' ' dR 2 

(3.36) 

where the total angular momentum J = P v is fixed. 
This is equivalent to the system 



dH 1 



elf 



-{r,p r = 0,p p ) = = 



Q2 H eS 



-{r,p r = 0,p v ) (3.37) 



or or 1 
with p v fixed, and this can be simplified to 

rWo(r 2 o + Jo) - ±r A' Q jl + 6A j 2 = 0, (3.38a) 



^o(r 2 +j 2 )-2A OJo 2 = 0. 



(3.38b) 



Here A' = dA(r)/dr, j = p v , and a subscript refers to 
quantities evaluated at the ISCO. The angular momen- 
tum can be solved for explicitly, leaving a single equation 
that must be solved numerically to determine the ISCO 
radius: 



r AoA% - 2roK) 2 + ^ A' = 0. 



(3.39) 



The angular orbital frequency of the ISCO is then found 
from Eq. (|3.27bD with p r = 0, 



dt 



'/'o-"o fl 



(3.40) 



For reference we also note that the EOB "horizon" is 
determined by solving Air) = . wh ile the EOB "light- 
ring" is found from the roots of jl2l| 



d (Air) 
dr 



0. 



(3.41) 



Logarithmic form of A(r) 



Building on previous works 122| - |l24 |. Barausse and 
Buonanno [94] have recently developed a new EOB 
Hamiltonian valid for spinning binaries. Their new 
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Hamiltonian has several interesting and desirable prop- 
erties. In particular, in the test-mass limit it reproduces 
the dynamics of the Mathisson-Papapetrou-Dixon equa- 
tions for a spinning point-particle in the Kerr spacetime 
(incorpor ating spin-orbit interactions to all PN orders) 
12iJl3lj |: and its PN expansion (for any mass ratio) re- 



produces the 3PN point-particle Hamiltonian, as well as 
the leading-order PN spin-spin interaction and the spin- 
orbit interaction up to 2.5 PN order. 

The details of this improved spinning EOB Hamilto- 
nian are quite complicated, but in the nonspinning limit 
the real and effective Hamiltonians match the forms given 
in the previous section. However, Barausse and Buo- 
nanno [94| have employed a different form for the func- 
tions A(r) and D(r) that appear in the effective metric. 
Rather than using Pade resummations of these functions, 
they introduce a logarithmic dependence which improves 
the behav ior o f A(r) in the spinning case. More specifi- 
cally, [9J, 1 1 32| found that when spins were present, the 
4PN and 5PN Pade versions of A(r) contain poles. Also 
the Pade resummation of A(r) did not always guarantee 
the existence of an ISCO in the spinning case, and when 
it did the ISCO did not vary monotonically with the spin 
magnitude. In the nonspinning limit, their new form for 
A(r) reduces to 

A l ° s (r) = (1 - vKy 2 [l - 2«(1 - T)K)} 
x [1 + log(l + Am + A 2 u 2 + A 3 u 3 + A 4 u% (3.42) 

where u = M/r, log refers to the natural loga- 
rithm, and the coefficients Ao through A4 are given in 
Eqs. (5.77)-(5.81) of [H| (with a = 0) and are functions 
of r\ and K. The function K = K(rf) parametrizes 4PN 
(and higher-order) corrections in Eq. (|3.42l) . It is given 
by [Eq. (6.11) of 0] 



K(ri) =K Q (l-4r 1 ) 2 +4(1 -277)77, 



(3.43) 



where the constant Kq = 1.4467 is chosen such that 
the resulting A l ° s (r) exactly reproduces the conservative 
GSF corrections to the ISCO computed by Barack and 
Sago Q. 



G. Shanks transformation 

A final method that we will consider is the "Shanks 
transformation," a nonlinear series acceleration method 
that can sometimes incre ase t he convergence rate of a 
sequence of partial sums [l 33| . This technique was in- 
troduced in the context of determining the ISCO in [o0| . 
The Shanks transformation relies on the approximation 
that the nth term in a converging sequence of partial 
sums Q n is related to the n — > 00 term Q by 



Q n = Q + ae r ' 



(3.44) 



with |e| < 1. By writing out equations for three succes- 
sive terms in the sequence (Q„_i, Q n , Q n +i) one can 



solve for the parameters Q, a, and e. Then for any ISCO 
quantity (e.g., the frequency, radius, or coefficient Cq de- 
fined below) with known values at 1PN, 2PN, and 3PN 
orders, we can define the Shanks transformation of that 
quantity by 



r)3PN,QlPN _ (f)2PN\2 
■\S _ V isco ^isco Wisco I 



S = 

^lSCO 



Q 



3PN 
isco 



)2PN 1 HlPN ' 
•> isco ' isco 



(3.45) 



This transformation will be applied to some of the ISCO 
methods discussed earlier in this section. 



IV. RESULTS 

For each of the methods reviewed in Sec. IIII| I have 
numerically computed the dimensionless angular orbital 
frequency Mf2 PN (?7) of the ISCO as a function of 77, and 
compared it with the renormalized Barack-Sago result 
[Eq. (|2.19p j. Specifically, I compute the analog of the 
coefficient Cq 11 



1.251 [Eq. (EH])] via 



PN 



lim — 

»7-KI 77 



n FN (77) 

iscoV // 



chw 



(4.1) 



where the limit is taken by evaluating at some sufficiently 
small value of r), and J7 PN is different for each method. 
The fractional error from the exact Barack-Sago result is 
also computed, 



„PN 



1. 



(4.2) 



Several of the methods discussed do not reproduce the 
standard Schwarzschild test-mass ISCO; these methods 
are ignored when computing Cq N . The remaining meth- 
ods for computing the ISCO are abbreviated as follows: 

1. C 2PN, C 3PN, C 4PN: the ga uge-invariant sta- 
bility condition from Sec. IIIIB1I at 2PN and 3PN 
orders. The pseudo-4PN version (Co4PN, not 
present in Table QJ fits a 4PN term to the exact 
BS result [see Eqs. (JB3J and flO]) below]. 

2. BfclPN, E h 2PN, E h 3 PN: th e hybrid energy- 
function method in Sec. IIIID H at each PN order. 

3. KWW-1PN, KWW-2PN, KWW-3PN: the Kidder- 
Will- Wiseman hybrid equations-of-motion ap- 
proach at each PN order (Sec. IIII D"2|) . 



4. HH-1PN, HH-2PN, 
Hamiltonian method 
PN order. 



HH-3PN : the hybrid- 
of Sec. IIII D 31 at each 



5. e2PN-P, e3PN-P: the e-method of Sec. HHE] using 
the 2PN and 3PN order Pade resummation of e(x). 

6. jlPN-P, j2PN-P, j3PN-P: the j-method of Sec. EH 
using the Pade resummation of j 2 (x) at each PN 
order. 
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7. A2PN-T, A3PN-T: the Taylor expansion of A(r) 
[Eq. ([531]) ] at 2PN and 3PN orders. 

8. A4PN-Ta, A4PN-T B : the 4PN Taylor expansion 
of A(r), with the two choices of the pseudo-4PN 
coefficient a 5 suggested in 

Choice A: A = 25.375, Ax = 0, and (4.3a) 
Choice B: A = -7.3, X 1 = 95.6. (4.3b) 

9. A2PN-P, A3PN-P, A4PN-P A , A4PN-P B , A4PN- 
Pc- analogous to items 7 and 8 above, but using 
the Pade approximants of A(r) listed in Eqs. (|3.35l) . 
A4PN-Pc7 uses a fit for the pseudo-4PN parame- 
ter 05 that exactly reproduces the BS conservative 
GSF ISCO shift [see Eq.flUE]) below]. 

10. AlogBB: denotes the logarithmic form of A(r) [94| 
given in Eq. (|3.42|) ; it is not listed in TableUJbecause 
it exactly reproduces the BS value by construction. 

11. HH-S, E h -S, KWW-S, and j-P-S all denote the 
Shanks transformation applied to the hybrid- 
Hamiltonian, hybrid energy-function, Kidder- Will- 
Wiseman, and j-methods, using the values for c™ 
for these methods at 1PN, 2PN, and 3PN orders. 

12. E1PN, E2PN, E3PN: uses the standard PN 
circular-orbit energy in Eq. (|3.1j) at 1PN, 2PN, and 
3PN orders to compute the ISCO. E-S denotes the 
Shanks transformation applied to the PN circular- 
orbit energy using the values from E1PN, E2PN, 
and E3PN. 

The resulting values for Cq and A cn are listed in Table 
HJ Figure [1] illustrates how some of the better-performing 
PN methods deviate from the exact BS value c™ 1 as a 
function of r\. Figure [5] shows the ISCO frequency for 
large values of r\ for several methods. 

Table [TT] lists the ISCO frequency in the equal-mass 
case for several of the methods discussed in Sec. lIIIi along 
with their fractional errors from the QCID results of [86| 
[rjg c c ID (l/4) « 0.1 2] 13 This value was chosen because 
(to my knowledge) Ref. [86[ seems to be the most re- 
cent and precise study of the ISCO using QCID calcu- 
lations. However it is not at all clear if this value ac- 
curately represents the "true" ISCO in the equal-mass 
case. This is especially true in light of the assumption of 
spatial conformal-flatness used in (86[; in a PN-context 
the spatial-metric is known to be conformally-flat only 
to 1PN order. One should thus interpret the PN com- 
parisons in Table [TT] with caution and with the under- 
standing that the exact value of the equal-mass ISCO is 




16 The values for the equal-mass ISCO in [8(1 vary from 0.121 to 
0.124 depending on the choice of method or boundary condition 
(the average is 0.122; see Table II of 86]). The comparisons in 
the second column of Table[ll]here drop the third uncertain digit. 



FIG. 1. (color online). The difference between various PN 
methods and the renormalized Barack-Sago conservative cor- 
rection to the ISCO frequency c™. We show only the meth- 
ods from Table [I] that have reasonable agreement (within 
~ 50%) with the exact value for c r n cn [Eq. fl20"|) ]. 



not accurately known (unlike the case of the BS conser- 
vative ISCO shift [l|, [3j). Nonetheless, I believe that the 
ISCO value quoted above from [86[ represents our cur- 
rent best-guess, so I will use that value in the remainder 
of this paper. 

Shortly after this article was accepted for publication, 
I became aware of an analy sis o f the ISCO using the 
"skeleton" approximation of [l34j . a truncation of Ein- 
stein's equations that assumes conformal flatness and 
drops some gravitational-field energy terms. The result- 
ing circular-orbit energy function that is derived from 
this approximation is computed to 10PN order; it agrees 
with the test-particle limit to all PN orders, but only 
agrees with the standard PN approximation to 1PN or- 
der for finite-^. The equal-mass ISCO frequency com- 
puted from this lOPN-order energy function (see Sec. VI 
of [H) is 0.0544, differing considerably from the 0.122 
value of [g6j]. I have also computed a hybrid version of 
this energy function (along the lines of Sec. IIII D the 
resulting value for c£ N (at the 10PN level) is -2.86, 
significantly different from the true value. Because these 
energy functions (and their lower PN-order variants) are 
based on a truncation of Einstein's equations and do not 
perform better than the top several approaches in Tables 
[TJand[lTJ I do not consider them further here. 



V. DISCUSSION 

From the values listed in Tables [TJ and [IT] we now make 
the following observations: 

1. Nearly all possible methods for computing finite- 
mass ratio corrections to the ISCO in the PN frame- 
work were considered, and these methods gener- 
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FIG. 2. (color online). ISCO orbital angular frequency as a function of the reduced mass ratio r\. The solid black line in each 
plot is the renormalized Barack-Sago result (extended to large rf). The abbreviations for each of the methods used are explained 
in the text. All EOB methods that use Pade approximants for A(r) are shown in the upper-left, and those that use Taylor 
expansions for A(r) (as well as the logarithmic approach of (94[) are in the upper-right. The lower-left plot shows the ISCO 
computed from different PN orders of the gauge- invariant ISCO condition in Sec. MI B II The lower-right plot shows results for 
the e-method, ji-method, and the 3PN hybrid energy function. The arrow near the point (0.25,0.12) indicates the equal-mass 
ISCO from quasicircular initial data calculations in [8(|. 



ally fall into two categories: nonresummcd and re- 
summed approaches. For the purpose of determin- 
ing the conservative ISCO shift for very small rj, 
all but one of the nonresummed approaches is use- 
less for computing c^ N since they generally do not 
reproduce the exact Schwarzschild ISCO. 

2. All of the methods discussed have appeared previ- 
ously in the literature, although not all were pre- 
viously investigated at 3PN order. In particular, 
note that the Kidder- Will- Wiseman j9l[ hybrid ap- 
proach (which originally motivated the develop- 
ment of resummation methods) was here extended 
to 3PN order. In contrast to the 2PN order re- 
sults reported in [9l|, at 3PN order the conserva- 
tive ISCO shift at least has the correct sign. How- 
ever, the fact that the 1PN version (KWW-1PN) 
makes the most accurate KWW prediction for c™ 



and the 3PN version (KWW-3PN) the least accu- 
rate, further suggests [H, HI that the KWW hy- 
brid approach is not a well-behaved resummation 
method. This pattern also occurs for the hybrid- 
Hamiltonian (HH) method (Sec. IIII D 3[) and the 
e-method (Sec. IIIIEf . suggesting that they too are 
not preferred approaches. This is in contrast with 
the remaining methods listed in Tables U and UH 
which share the property that the higher PN iter- 
ation of a given method produces a value closer to 
the exact result. 

3. The method that produces the best agreement with 
the exact result (~ 10% error) is the EOB method 
in which a pseudo-4PN parameter 05(77) is intro- 
duced and its value is adjusted to NR simulations 
in [Hi]]. In p articular, only one of the two suggested 
choices [5J] for 05 [choice A in Eq. (|4.3[) ] gives good 
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TABLE I. Tabulation of the conservative GSF correction to 
the ISCO for several of the PN methods presented in the text. 
The first column lists the PN method used (descriptions of the 
abbreviations are given in Sec. IIV[) . The second column lists 
the coefficient Cn N defined in Eq. (|4.1I) . This was computed 
for r\ — 1CP 6 ; for smaller r/ the values are unchanged to four 
significant figures. These values are compared with the exact 
(renormalized) Barack-Sago result, cff n = 1.251, in the third 
column, which lists the fractional error [Eq. (|4.2[) ]. The table 
is sorted by the absolute value of the fractional error in the 
third column (most accurate method first). 



Method 






A4PN-Pa 


1.132 


— 0955 


A4PN-Ta 


1.132 


-0.0955 


C 3PN 


1.435 


0.1467 


e2PN-P 


1.036 


—0.1717 


KWW-1PN 


1.592 


0.2726 


A3PN-P 


0.9067 


-0.2754 


A3PN-T 


0.9067 


—0.2754 


A4PN-P B 


0.8419 


—0.3272 


A4PN-T S 


0.8419 


—0.3272 


j3PN-P 


1.711 


0.3671 


J2PN-P 


0.6146 


—0 ^088 


KWW-S 


0.5610 


-0.5515 


C* 2PN 


0.5833 


-0.5338 


E h 3PN 


0.4705 


-0.6240 


e3PN-P 


2.178 


0.7409 


A2PN-P 


0.2794 


-0.7767 


A2PN-T 


0.2794 


-0.7767 


E h 2PN 


0.0902 


-0.9279 


E h lPN 


-0.014 73 


-1.011 


E h -S 


-0.054 71 


-1.044 


HH-S 


-0.1486 


-1.119 


jlPN-P 


-0.1667 


-1.133 


KWW-2PN 


-1.542 


-2.232 


j-P-S 


-2.104 


-2.682 


KWW-3PN 


4.851 


2.877 


HH-1PN 


6.062 


3.844 


HH-2PN 


-12.75 


-11.19 


HH-3PN 


25.42 


19.32 



agreement. Choice B gives an error that is 3 times 
worse. It is especially interesting that the fits to 
the NR simulations — which are done in the q ~ 1 
limit — have in some sense "preselected" a value for 
05 that is closest to reproducing the exact result of 
a<j< 1 calculation. 

4. It is interesting to note that if we neglect the meth- 
ods that involve some sort of fitting to numerical 
results, then the best EOB approach (A3PN-P) is 
not the most accurate method. Rather, in both the 
extreme-mass ratio (Table HJ and equal-mass (Ta- 
ble[ITJ) cases, two distinct nonresummed approaches 
based on the ordinary 3PN equations of motion are 
among the most accurate approaches. 

5. However, note also that in both the extreme-mass 
ratio and equal-mass cases, the error associated 
with the A3PN-P method is nearly the same (~ 



TABLE II. ISCO frequency for equal-mass binaries for se- 
lected methods presented in the text. The first column lists 
the PN method used (descriptions of the abbreviations are 
given in Sec. IIV[) . The second column lists the ISCO angu- 
lar orbital frequency MQi aco (rj = 1/4). The third column 
lists the fractional error from the approximate QCID value 
Q ?sco D ~ 0.12 reported in f8f|. The table is sorted by the ab- 
solute value of the fractional error listed in the third column. 



Method 


a 6 lSCO 


a QCID 


j3PN-P 


0.1207 


0.0061 


E-S 


0.1285 


0.071 


E3PN 


0.1287 


0.073 


e3PN-P 


0.1340 


0.12 


A4PN-P C 


0.1036 


-0.14 


E2PN 


0.1371 


0.14 


A4PN-P.4 


0.1004 


-0.16 


A4PN-P B 


0.098 07 


-0.18 


AlogBB 


0.089 99 


-0.25 


e2PN-P 


0.088 50 


-0.26 


A3PN-P 


0.088 22 


-0.26 


C 4PN 


0.1567 


0.31 


C 2PN 


0.0809 


-0.33 


j2PN-P 


0.079 80 


-0.33 


E h 3PN 


0.076 98 


-0.36 


A2PN-T 


0.073 40 


-0.39 


A2PN-P 


0.073 12 


-0.39 


E h 2PN 


0.069 59 


-0.42 


£UPN 


0.067 79 


-0.44 


E h -S 


0.067 21 


-0.44 


jlPN-P 


0.065 30 


-0.46 


j-P-S 


0.057 35 


-0.52 


E1PN 


0.5224 


3.4 



27%), and arises from a single, distinct method. 
Introducing a pseudo-4PN term and calibrating to 
NR (A4PN-P A ) or to the BS results (A4PN-P C ) 
further reduces the errors in both mass-ratio lim- 
its. 

6. In the equal-mass case (where the nonresummed 
PN series has good convergence properties), the 
ISCO computed from the 3PN circular-orbit en- 
ergy shows remarkable agreement (~ 7%) with 
the QCID result from 0— better than any EOB 
method. 

7. The j-method at 3PN order (j3PN-P) produces 
nearly exact agreement with the equal-mass QCID 
ISCO. This is possibly coincidental, and partly due 
to the truncation of the "exact" QCID result to 
two digits. However, we also note that j3PN-P does 
moderately well at reproducing the BS ISCO shift 
(~ 37% error), and the different PN iterations of 
the j-method (1PN, 2PN, 3PN) show successive im- 
provement at each PN order in both the equal-mass 
and extreme-mass-ratio cases (in contrast with the 
e- method; see Tables Q] and [TTJ) . 

8. Surprisingly, the method which best reproduces the 
Barack-Sago conservative GSF ISCO corrections 
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(without relying on any fits with NR or GSF cal- 
culations) is the gauge-invariant ISCO condition 
Co3PN of (62[. This is the only nonresummed 
approach in Table |U It is an especially interest- 
ing method because it both reproduces the exact 
test-mass ISCO and matches the conservative GSF 
ISCO corrections with good accuracy (~ 15%) — 
nearly twice the accuracy of the best 3PN EOB 
approach. On the other hand, the Co3PN method 
does not produce an ISCO for -q > 0.183, while the 
A3PN-P method yields a well-defined ISCO for any 
V < 1/4. 

9. Notice also that in Table HI the Taylor and Pade 
forms of the EOB potential give identical results 
at a given PN order. This arises from the fact [dis- 
cussed after Eq. p.35g| )] that the 77-expansion of the 
Pade approximants of A(r) reduces to the Taylor 
expansion given in Eq. (|3.3ip . 

10. Reg arding the logarithmic form of A(r) proposed in 
[94| . note that the AlogBB ISCO frequency main- 
tains a nearly constant slope for all 77, closely fol- 
lowing the large-77 extrapolation of the BS result to 
which it is calibrated (Fig. [2j . 

1 1 . Regarding the Shanks transformation: while it was 
only possible to apply it for methods with ISCO 
quantities defined at 1PN, 2PN, and 3PN orders, 
for the methods investigated in Tables U and UU it 
generally did not yield an improvement in accu- 
racy. The exception is for the standard PN-energy 
function method E3PN, which saw a very slight im- 
provement in the accuracy of the equal-mass ISCO. 

12. Also note the much larger spread in the error val- 
ues listed in Table U versus those in Table |TTJ This 
is likely a reflection of the well-known poor con- 
vergence of the PN series in the small-77 limit and 
the relatively good convergence in the equal-mass 
limit. 

Let us now elaborate on point|8]above. It is rather curi- 
ous that the gauge-invariant ISCO condition of Blanchet 
and Iyer [62| [Eq. (|3.7|l ] not only exactly reproduces the 
test-mass ISCO, Xi sco = 1/6, but also provides very close 
agreement with the BS result for Cjf n . Damour [71j also 
finds a value for c™ equal to that produced by Co3PN; 
but his result is derived by performing a PN expansion 
of quantities that appear in the EOB effective metric, 
so the agreement with the Schwarzschild ISCO in his ap- 
proach is by construction [see his Eq. (5.41) and the asso- 
ciated discussion]. In Blanchet and Iyer's [62| derivation, 
the agreement with the test-mass ISCO was not enforced 
(and is thus surprising since standard PN calculations 
typically do not exactly reproduce strong-field results). 



As for the fact that the value c, 



C 3PN 



1.435 [from 

Co3PN or Damour's [7l| Eq. (5.41)] agrees more closely 
with the BS result than the A3PN-P value, Damour [71| 
suggests that this is an accident arising from the failure of 



certain terms in the PN expansion of EOB quantities to 
cancel with higher-order terms. However, it is not clear 
how (or if) this argument translates to an "explanation" 
of the value of c^° 3PN ~ 1.435 when it is derived via the 
standard PN approach in Blanchet and Iyer [62| . 

To further explore the possibility that the behavior of 
the gauge-invariant ISCO condition Co3PN might be ac- 
cidental, I have extended the calculation of Blanchet and 
Iyer [62| to the case of nonprecessing, spinning binaries. 
The details are presented in (95|. The derivation follows 
that in [62l |91| , except that I als o include the spin-orbit 
terms at 1.5PN and 2.5PN orders [IMIIqI, and the spin- 
spin and quadrupole-monopole terms at 2PN order [l35j j . 
The result is a condition for the ISCO which generalizes 
Eq. (JX71), 



s _ M 2 
Co = — o-Cq — 1 



T 3 

.1 Q 



6iEn 



3/2 



14 



M 2 



5m 



U77-3 



M 2 ) 



5/2 



'"'(22 + 32,) -^§(18 + 15,) 



M 2 



/397 123 ,\ 2 



(5.1) 



where Sf = I ■ (SJ + E^j = Mi ■ {S^/m 2 - Sf/mi), 
= M£-[(l+m 2 /m 1 )S^ + {l+m 1 /m 2 )S 2 z ),eis the unit 
vector in the direction of the Newtonian orbital angular 
momentum, and 5m — mi — m?P^l 

If the larger BH has spin \S%\ — X2 m 2 ( an d Si — 0), 
the above condition reduces, in the test-particle limit, to 



Co = 1 - 6x0 + 8 X ^o /2 



3(x! 



A2^2 



4x c 2 4 /2 



0{x%). 
(5.2) 

A similar criterion for the exact (asymptotic) Kerr ISCO 
frequency parameter X = (m 2 ^°^ ) 2 ^ 3 can be derived 
from the minimum of th e orbital energy of a particle in 
the Kerr spacetime [l36j ] . This criterion, when expanded 
for small BH spin, takes the form [95j] 



CT 1 ' 1 = 1 - 6A + X2 (8A 3/2 - 4A7T) 

+ xl (-3X 2 + 8A 3 - 10A 4 /3) + 0(xD- (5.3) 

Note that no PN expansion was used to derive this ex- 
pression; it is exact to 0(x|) m t ne BH spin. 

Except for the 3PN and 4PN order spin terms (whose 
forms in the PN equations of motion are not currently 
known), Eqs. (|5.2p and (|5.3p agree exactly when we 
identify the Kerr spin parameter X2 with xl (note that 



'5/2x 



17 The spin angular-momentum vecto rs Sf used ab ove refer to the 
constant-magnitude spin variables [104, 106, 13EJ. However, the 
same test-mass limit [Eq. 15.20 ] is found if the spin variables with 
varying magnitudes 105, IO3] are used (see I95|l '). 
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xq — > Xq in the test- mass limit). The fact that the exten- 
sion of the gauge-invariant ISCO condition Co3PN [62] 
to spinning BHs reproduces the exact Kerr ISCO (to the 
expected order in \2) suggests (but does not prove) that 
the close agreement of Co3PN with the conservative GSF 
ISCO shift is not accidental. It remains to be seen if the 
predictions of Eq. (|5.1[) will also produce good agreement 
with the conservative GSF ISCO shift when it is even- 
tually calculated for Kerr BHs. In addition to this good 
agreement with the Kerr ISCO, condition (I5.1[) can also 
be shown to exactly reproduce the fully-relativistic con- 
servative shift in the ISCO due to the spin of a small 
test- mass. These issues are discussed further in 19511 . 



VI. JOINT CONSTRAINTS BETWEEN 
NUMERICAL RELATIVITY, 
POST-NEWTONIAN THEORY, AND 
SELF-FORCE CALCULATIONS 

The primary motivation for GSF calculations is to pro- 
duce accurate waveform models for EMRIs. However, 
as we explore below, GSF calculations can also improve 
our understanding of comparable-mass (1/10 g ;$ 1) 
waveforms. For example, GSF calculations can help cal- 
ibrate EOB waveforms (Sec. |VTA1 see also (26l. FFlL 19^). 
phenomenological inspiral-merger-ringdown (IMR) wave- 
forms ([Sec. lVlBI . and high-order terms in PN quantities 
like the orbital energy (Sec. IVTC]) . In Sec. IVIDl I also 
discuss how numerical computations of sequences of qua- 
sicircular initial data, combined with full NR calculations 
of the energy flux for a few orbits, could be used to esti- 
mate high-PN-order terms in the inspiral phasing. 

A. Fitting pseudo-4PN parameters 

In addition to the two choices for as in Eq. (|4.3I) . we 
can also attempt to fix the value of 05 such that the 
Barack-Sago conservative GSF ISCO shift is exactly re- 
produced. To do this we assum^H that Ai = (since 
the GSF calculations are only to leading-order in 77), and 
adjust A (holding r\ fixed at 10~ 7 ) until the resulting 
matches the exact result. This is done for both the 
4PN Taylor and 4PN Pade expansions of A(r) (denoted 
A4PN-T C and A4PN-P C here and in Figs. Q] and HJ}. The 
resulting value for Ao is 

Aq* = 38.84 (6.1) 

for both A4PN-T C and A4PN-P c 0f| Note that this value 
is not wildly different from Aq = 25.375, one of the values 



We also ignore the possibility of any log terms that might be 
present at 4PN order; see [26l, |71|| for further discussion. 
One can also attempt to leave undetermined the 3PN coefficient 
a 4 [Eq. JOS}] and "fit" it to cff n . The result is a 4 « 28.96r?, an 
~ 55% error with the exact result. 



found by fitting to the NR simulations in Ref. [54j [see 
Eq. (|4.3|) ]. Damour [7l[ finds an equivalent constraint 
[see his Eq. (4.42)], although he considers a 5PN exten- 
sion of the EOB potential A(r) that supplements the 3PN 
Taylor expansion by 8 A = rj(a^u 5 + a^u 6 ) [50j, where 
u = M/r and a§ and a§ are constants. 

The above fit for the pseudo-4PN parameter as = A^yy 
does not necessarily contain any new physical informa- 
tion about the 4PN dynamics; instead it essentially "re- 
sums" all of the higher-order PN corrections that con- 
tribute to the conservative GSF ISCO shift and groups 
them into a single 4PN term. This value for 05 would not 
exactly coincide with some future calculation of the 4PN 
expansion of A(r); however, if the 5PN and higher con- 
tributions to A(r) are small, it is possible that the value 
calculated above might not be too far from the true 4PN 
result. 

To test this notion, compare the 0(rf) terms in the 4PN 
Taylor expansion of A(r) [Eq. (|3.31[) ]: 2r//r 3 , a,i/r A , and 
a§/r° . Numerically evaluating these terms at r = 6 [and 
using Eq. dHUJ for a 5 ], we get 0.0093??, O.OUrj, and 0.005?? 
(respectively). The fact that the 4PN term is smallest 
even near the ISCO suggests (but does not prove) that 
the 5PN term might be small in comparison to the 4PN 
one. See Damour [Tljj for an alternative argument. 

Regardless of its agreement with respect to the true 
value, the above choice [Eq. (|6.1[> ] for the pseudo-4PN 
term in A(r) can be said to accurately reproduce an im- 
portant strong-field feature of the conservative dynamics. 
It might therefore be useful to fix Ao to the value 38.84 in 
future EOB/NR comparison studies. Fixing the pseudo- 
4PN term in this way could also be useful for studies 
that use the EOB formalism to model waveforms from 
extreme- and intermediate- mass-ratio inspirals [64| . 

If one also had a highly accurate determination of the 
ISCO frequency in the equal-mass case (or any moderate 
mass-ratio for that matter), one could also attempt to 
constrain higher-order parameters in the 05(77) coefficient 
[Eq. (|3.34[) ]. For example, Ref. [86[ gives values for the 
nonspinning, equal-mass ISCO in the range Mfi(l/4) = 
0.121 to 0.124 depending on the method and boundary 
condition used. If we fix the value of Ao to that found 
in Eq. (|6.ip . then we find that the above range for the 
ISCO implies that Ai must lie in the range 

Ai <= [488.9,665.5]. (6.2) 

For MQ(l/4) = 0.122 we find Ai = 541.3. 

In addition to using the BS conservative ISCO shift, 
Damour [Til ] also used comparisons with full NR simu- 
lations to provide constraints on the parameters a§ and 
a® discussed above. Rather than using full NR simula- 
tions (which include the dissipative dynamics), sequences 
of BH quasicircular initial data (combined with the con- 
servative GSF ISCO correction) could also be used to 
constrain these parameters [as was done with (Ao,Ai) 
above]. This has the advantage of using strictly "con- 
servative" NR calculations to constrain the conservative 
dynamics encapsulated in A(r). 
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We can similarly try adding a_pseudo-4PN term to 
the gauge-invariant Blanchet-Iyer |62j stability condition 
discussed in Sec. IIIIB 11 



C = 1 - 6a;o + Ur/xl 



397 



123 



r] - 14r/ 2 ) Xq + c 4PN t]Xq, 



(6.3) 



where we again ignore possible logarithmic corrections 
(computable from the results of [73|]) and assume that 
the 4PN terms do not modify the stability condition at 
0{rf) (allowing the exact Schwarzschild test- mass ISCO 
to be reproduced). Precise agreement with the Barack- 
Sago result is enforced by tuning C4pn to the valuJ^l 



c 4PN ps -158.64. 



(6.4) 



since the smallest mass ratio considered in [l37j j was 0.25 
(r/ = 0.16); their fits could only be expected to work for 
mass ratios larger than this. Future phenomenological 
IMR templates could consider fixing the value of y( 10 ) 
to 6- 3 / 2 d£ n ps 0.085 14; this might help to provide bet- 
ter template matches with NR simulations at small mass 
ratios. In principle the other higher-order r] terms in 
Eq. (JHUD could also be fixed via any of the PN ISCO 
methods discussed here (see, e.g., Fig. [2]) or by QCID 
ISCO calculations. However, it is not clear if doing so 
will necessarily produce a template family that can bet- 
ter match NR waveforms. 



C. Constraints on the 4PN and 5PN circular-orbit 

energy 



The resulting ISCO frequency is denoted Co4PN in Ta- 
ble [XT] and Figs. Q] and [5] Evaluating the 0(ij) terms in 
Eq. (O) at x ps 1/6, 



Urjx 2 ps 0.397?, 



397 



123 

7 
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c 4PN r)x 



0.57r?, 
-0.127?, 



(6.5a) 
(6.5b) 
(6.5c) 



we see a similar sequence as in the EOB case with the 
pseudo-4PN term being smallest. 



B. Constraints on phenomenological 
inspiral-merger-ringdown templates 

In the phenomenological inspiral-merger-ringdown 
template family developed by Ajith et al. [1371 ] . a 
frequency-domain template h(f) = ^4(/)e~ 1 *^ is de- 
fined in te rms o f a phase $(/) and an amplitude A(f) [see 
Eq. (1) of [137]]. The amplitude function A(f) is written 
as a piecewise function that transitions from an "inspi- 
ral" to "merger" at a frequency /i , and to a "ringdown" 
a t f%- The frequency fi = Qi/tt between the inspiral and 
merger phase is defined such that in the i] — > limit it 
reduces to the test-mass ISCO. When 77 is nonzero and 
for initially nonspinning BHs, this tr ansit ion frequency is 
given by [see Eq. (2) and Table I of [US]] 



MOi = 6~ 3/2 



(10) 



r?V 20) 



7? 3 y (30) i 



(6.6) 



where y^ = 0.6437, = -0.058 22, and y^ ~> = 

-7.092. Note that 6 3 /V 10) ps 9.46 and differs from cjf 1 
by a factor of about 7.6 (although it at least has the 
correct sign). This lack of agreement is not surprising 



The PN energy for circular orbits [Eq. (|3.1I) ] contains 
newly computed terms [lU at 4PN and 5PN orders. How- 
ever, only the test-mass-limit and logarithmic pieces of 
these terms are known. The remaining uncertainty is 
parametrized by the polynomials 

e 4 (ry) - ef + e 4 % + eft? + • • • efrf, (6.7a) 
e 5 (v) = 4 0) + 4 1} V + 4V + • • • 4V , (6.7b) 

where p and q are integers (probably equal to 4 and 5). 
Computing these functions will require the completion of 
the PN iteration scheme at the 4PN and 5PN levels — a 
daunting task. Here I attempt to partly constrain these 
polynomials by using numeric calculations of the ISCO 
from the GSF approach or quasicircular initial data cal- 
culations. 

Numerical relativity calculations of the ISCO fre- 
quency for equal-mass binaries come from examining 
QCID sequences. The most recent results [86| indicate 
Mf2i SCO (l/4) pa 0.122. This single value provides a possi- 
ble crude constraint on PN parameters appearing in the 
circular orbit energy function through the ISCO condi- 
tion d/dx[E PN (fl = fi isco (l/4),ry = 1/4)] = 0, where 
S PN = E PN /(r]M) is given by Eq. j3l[])E3 We can test 
this approach by assuming that the 3PN coefficient in 
E FN is unknown and has the form 



675 

where the known value for 63(77) is 



e3(f/0 



34445 205 



96 



576 



155 
"96" ? 



35 

7 

5184 



(6.8) 



(6.9) 



If we perform a similar procedure leaving the 0{rf) 3PN coeffi- 
cient in Eq. 113.711 undetermined, we find a value C3PN ~ 96.2, 
which is within 22% of the exact result. 



Keep in mind that it is not clear to what degree the QCID results 
represent the "true" ISCO embodied by the purely conservative 
dynamics of the full Einstein equations. This is especially true 
in light of the conformal-flatness assumption that is used in 86] . 
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Ignoring the 4PN and higher order terms and using the 
above ISCO criterion with the equal-mass value from [86[ 
yields the constraint 

e 3 (l/4) » 34.4, (6.10) 

which agrees with the exact value of 38.3 to 10%. If we 
apply this procedure to the 4PN and 5PN corrections, 
we get the constraints 

64(1/4) 147.3 at 4PN order, and (6.11a) 

e 4 (l/4) + 0.2952e 5 (l/4) ps 189.6 at 5PN order. 

(6.11b) 

The reliability of these constraints is not completely 
clear; it depends on the accuracy and precision of the 
numerical simulation (including the systematic effects 
alluded to above), as well as on the contributions of 
all the higher-order PN terms. For example, suppose 
that the ISCO computed from [§6| has a fractional er- 
ror of 0.002/0.122 ss 2%. Any PN corrections that one 
might hope to resolve should contribute at least 2% to 
the value of the ISCO. Using the test-mass limit as a 
gauge, the nth-order PN expansion of Eq. (|3.2[) pre- 
dicts fractional errors for the Schwarzschild ISCO of 
^nPN/g-3/2 _ 1 _ [7.0,0.82,0.27,0.12,0.055,0.028] for 
n = 1 to 6. This indicates, for example, that the 5PN 
and higher-order terms account for 12% of the ISCO, 
6PN and higher-order terms account for 5.5%, etc. This 
suggests that the simulations in [86[ should be precise 
enough to resolve 4PN and 5PN effects, although not 
necessarily with high accuracy or without contamination 
from higher-PN terms. 

We can attempt to also derive constraints on 64(0) 
and 65(0) using the BS conservative GSF corrections 
to the ISCO. In this case we compute the condition 
d/dx[E^^ lid (Q — Ofl^t Sf C 7 ?))] = 0, where we substitute 
Eq. (|2.19p for the frequency into a hybrid energy func- 
tion analogous to Eq. (pT22|) [but with the 4PN and 5PN 
terms in Eq. (|3.ip included]. The resulting ISCO condi- 
tion is then expanded to linear order in 77, yielding the 
constraints 

e 4 (0) w 429.1 at 4PN order and (6.12a) 
e 4 (0) + e 5 (0)/5 382.8 at 5PN order. (6.12b) 

To further constrain the functions 64(77) and 65(77) we 
make the following two observations: first, if we exam- 
ine the numerical values of the coefficients in the 3PN 
expansion of E PN [cf. Eauation (j3.ip ]. 

l+a;(-0.75-0.08377)+x 2 (-3.4+2.4r7-0.042r7 2 ) 

x 

+ x 3 (-ll. + 39.77 - I.677 2 - O.OO6877 3 ) + 0(x 4 ), (6.13) 

we see that the coefficients of each power of 77 tend to 
decrease in absolute value as the power of 77 increases. 
Second, since 77 is at most 0.25, the absolute value of 



the terms proportional to rf is further suppressed by a 
factor 0.25 p or smaller. This suggests that we can ap- 
proximately ignore some of the higher-order 77-terms in 
the expansions in Eqs. (|6.7[) . For example, working at the 

4PN level only, we can approximate 64(77) ps e^ ' + r\e^ . 
Then Eqs. (|6Tla|) and (|6.12ap imply 

e 4 0) « 429.1 and ef 1 ps -1127. (6.14) 

One could do something similar at the 5PN level, but 
as there are more unknown parameters than equations, 
one would need numerical values of the ISCO for more 
values of 77. These values should be computable via QCID 
calculations analogous to those in j86[, and could also 

help to constrain the higher-order coefficients ef\ 

etc. One might also suspect that the terms ef^ and 
could be constrained by current GSF calculations as was 
done with the redshift function u\ (y) in [73| . 

The estimates on 64(77) and 65(77) presented here are 
meant to illustrate techniques through which they could 
be constrained. Better constraints would require more 
accurate numerical simulations for several mass ratios. 
Nonetheless, if the exact 4PN and 5PN terms are even- 
tually computed, it would be interesting to compare their 
values with the above estimates. 



D. A suggested approach for numerically 
computing higher-order PN corrections to the 
gravitational-wave phasing 

The above subsections indicate that the recent GSF 
results for the conservative shift in the ISCO can better 
inform our knowledge of comparable-mass waveform tem- 
plates. However, as they are currently limited to the first- 
order self-force, GSF calculations are constrained to only 
provide information about the q <C 1 limit of these tem- 
plates. Numerical relativity provides exact comparable- 
mass waveforms, but becomes severely limited by com- 
putational costs for q < 1/10. These costs are especially 
severe if one desires many cycles before the merger (in the 
regime where our current 3PN waveforms are starting to 
lose phase coherence). Here, I propose an alternative in- 
spiral template-generation method based on calibrating 
higher-order PN terms with low-cost (or at least "lower- 
cost") NR simulations. 

Match-filtering-based detection and analysis methods 
are most sensitive to the phase of the gravitational wave- 
form. For quasicircular binaries, the phase of the Fourier 
transform of the GW signal ^f(f) is determined by two 
ingredients from PN theory, the orbital energy E{x) and 
the energy flux [actually the luminosity C gw (x) = —E] 
as a function of the orbital frequency MQ = x 3 / 2 , 

d 2 *(f) 9 dE/df 

-^ = - 2 *c^rr (6 - 15) 

where / = is the GW frequency. For nonspin- 

ning binaries, these two ingredients are currently known 
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to 3.5PN order. Computing the 4PN and higher-order 
terms will be very difficult, and it is not clear if they will 
be computed before the first detections. While the 3.5PN 
order terms are sufficient for detecting GWs, knowledge 
of higher-order phase corrections will allow the recovery 
of more signal-to-noise at later times into inspiral, allow- 
ing expensive NR simulations to focus on the cycles very 
close to merger. 

Rather than run standard NR evolutions to compute 
waveform templates for "smallish" (q ~ 1/10) mass ra- 
tios, I instead propose the following two-pronged strategy 
that involves using NR to calibrate undetermined param- 
eters in standard (non-EOB) waveform templates. The 
first part involves computing higher-order corrections to 
the orbital energy E(x). As shown above, the 4PN and 
5PN logarithmic pieces of this function have been re- 
cently computed [73|], but the nonlogarithmic terms are 
unknown. Using quasicircular initial data computations 
of the equal-mass ISCO and the GSF conservative ISCO 
shift, Sec. IVI CI set some additional constraints on the 
undetermined functions at 4PN and 5PN orders. My 
primary suggestion is to use several QCID calculations — 
at a variety of mass ratios — to set further constraints on 
the undetermined 4PN and 5PN pieces of E(x). While 
this could be accomplished by computing the ISCO for 
a variety of mass ratios, a better strategy might be to 
compute the energy for a variety of frequencies and mass 
ratios. This would essentially determine E(x) up to 4PN 
or 5PN order through numerical fits to the QCID results. 
This is analogous to the fitting procedure used in [73[ 
to determine high-PN-order terms in the gauge-invariant 
redshift function u\ . If run for small enough mass ratios 
(<Z ^ 1/100), QCID calculations could potentially provide 
additional points of comparison with GSF calculations. 

The second step is to "compute" the GW luminos- 
ity to higher PN order than 3.5PN. In this case we are 
partially helped by analytic BH perturbation theory cal- 
culations of the test-mass limit terms in the luminosity, 
which are currently known to 5.5PN order (see jl38j and 
references therein). It might also be possible to extend 
the program of Ref. [73] to the computation of the finite- 
mass-ratio logarithmic terms in the luminosity at 4PN 
and higher orders. To obtain the remaining finite-mass 
ratio nonlogarithmic terms, one fits these undetermined 
coefficients by comparing with the luminosity from full 
NR evolutions. These evolutions are very expensive if 
the inspiral starts at large separations or has small mass 
ratios. However, to fit higher-PN terms in the luminosity 
we do not necessarily need full evolutions of the entire in- 
spiral and merger. Rather, we could make do with small 
stretches of a simulation that consist of only a few orbits 
near a single frequency. This allows us to improve the 
fit to C gw (x,r]) by supplementing the currently available 
NR values of £ gw with a few discrete points at large sep- 
arations and/or small mass ratios. Although even these 
few-orbit evolutions might still be expensive, they could 
have a big payoff in providing a permanent calibration of 
the PN phasing [determined via Eq. (|6.15[) ]. 



In practice, there are several difficulties associated with 
the above scheme. One obvious issue concerns the accu- 
racy of current QCID calculations. Many (but not all) 
of these calculations assu me that the spatial metric is 
conformally-flat (see [l39j and references therein); this is 
known to be accurate only to 1PN order. While some 
higher-order PN effects are still implicit in these calcu- 
lations (which indeed show good agreement with 3PN 
calculations, cf. Table ITU . one would clearly like a calcu- 
lational scheme that is at least self-consistent to the order 
of the PN corrections that one is trying to compute (in 
this case 4PN). In computing the GW luminosity, there 
are also problems associated with performing NR simu- 
lations at large separations (and small mass ratios) : evo- 
lutions are slow for these orbits, a large computational 
grid is required, and one must be careful of contamina- 
tion from junk radiation and boundary reflections. One 
must also perform the calculations at separations large 
enough that the adiabatic approximation [upon which 
Eq. (I6.15[) relies] is valid. However, at large separations 
the higher-order PN terms that one is trying to fit also 
become increasingly small, making them potentially dif- 
ficult to resolve in a numerical simulation. 

To gauge the needed precision, we can consider the 
size of the 4PN terms of E{x) and £ gw in the test-mass 
limit. The fractional error in the size of the 4PN term, 
(E 4PN - E 3PN )/E Schw varies from ~ 0.0002 at x = 1/20 
to - 0.004 at x = 1/10 to - 0.03 at x = 1/6. For 
comparison, I estimate that the ISCO energy calculated 
in [86j has a precision of roughly 1% (see their Table II). 
This is just enough precision to resolve the 4PN term near 
the ISCO, but not enough to resolve it at much larger 
separations. Similarly, we can examine the PN expansion 
of the GW luminosity, which is known to 5 .5PN order in 
the test-mass limit [see, e.g., Eq. (174) of [l38j ]. In this 

varies 



/*3.5PN 

gw 



)/£ 5.5PN 



case the fractional error (£ g ™ 

from ~ 0.001 at x = 1/20 to ~ 0T02 at x~= 1/10 to 
~ 0.1 at x = 1/6. Current NR simulations can attain 
this level of p recision out to at least x ~ 1/10 (see, e.g., 
Fig. 2 of [38|). Future work will examine in more detail 
the feasibility of the scheme proposed here. 



VII. CONCLUSIONS 

The primary purpose of this study has been to com- 
pare the recent gravitational-self-force (GSF) calcula- 
tions of the conservative shift in the Schwarzschild 
ISCO to nearly all PN/EOB methods for computing the 
ISCO. The results, summarized in Table |H show that 
while EOB methods calibrated to NR simulations per- 
form best, uncalibrated EOB — as well as other resumma- 
tion approaches — do not perform better than the gauge- 
invariant ISCO condition of 



62]. This ISCO condition 



has the especially interesting property of exactly repro- 
ducing the Schwarzschild ISCO, even though it is derived 
using the standard PN equations of motion without any 
form of resummation (which is typically used to enforce 
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the test-mass limit). To investigate if this agreement is 
accidental, I have also generalized this gauge-invariant 
ISCO condition to spinning BHs, and showed that it re- 
produces the Kerr ISCO up to the expected order in the 
spin and PN expansion parameters. This approach also 
exactly reproduces the fully-relativistic conservative shift 
in the ISCO due to the spin of the test-mass (see [95[ for 
details). 

The various PN/EOB ISCO methods were also com- 
pared with the quasicircular initial data calculations of 
the equal-mass ISCO in [86| . In this case, a nonresummed 
PN method also performs better than EOB approaches. 
However, this is a different method (one based on the 
minimum of the 3PN orbital energy) than the Blanchet- 
Iyer [H3] ISCO condition. The 3PN-EOB approach has 
the advantage of being a single resummed method that 
can model the ISCO in both the comparable-mass and 
test-mass limits with comparable (albeit larger) errors in 
both cases. 

These results suggest that the standard PN equa- 
tions of motion somehow contain information about the 
strong-field conservative dynamics (at least in the test- 
mass limit). This is surprising since PN quantities often 
converge slowly in the rj — > limit. The gauge-invariant 
ISCO condition of [62| (and its generalization to spinning 
binaries) apparently does not suffer from this limitation. 
Subsequent work will examine predictions of the conser- 
vative ISCO shift in Kerr from EOB and PN approaches 
[95( 1 . These predictions can be compared with future GSF 
calculations in the Kerr spacetime. 

The « 28% error between the 3PN EOB prediction for 
the conservative GSF ISCO shift and the exact Barack- 
Sago result suggests that while the EOB formalism ex- 
actly encapsulates the test-particle limit of motion in 
Schwarzschild, it does not encapsulate small-deviations 
from the test-mass limit any better than standard PN 
approaches. While this is not necessarily unexpected, 
it suggests caution be used when attempting to model 
conservative GSF effects in EMRI waveforms using EOB 
methods (64|. Of course, the EOB formalism can be mod- 
ified via the introduction of unknown terms that can be 
calibrated to GSF calculations [26], [7l[ or to numerical 
relativity. However, for intermediate mass ratios there 
is no accurate numerical method to calibrate against, so 
it is useful to compare the performance of different ap- 



proaches in the absence of any calibration. 

This study has also explored how GSF calculations 
can further our knowledge of comparable-mass template 
waveforms. GSF calculations of the conservative ISCO 
shift can be used to calibrate parameters in EOB and 
phenomenological IMR waveforms, and, combined with 
quasicircular initial data calculations of the comparable- 
mass ISCO, can help constrain the undetermined func- 
tions in the 4PN and 5PN pieces of the PN orbital energy. 
Calculations of quasicircular initial data sequences at un- 
equal mass ratios are needed to further constrain these 
parameters and functions. 

A new method of calibrating 4PN (or higher) terms in 
the waveform phasing was also proposed. Rather than 
comparing with the full-NR waveform phase, this ap- 
proach suggests using two separate NR calculations to 
determine the phase at a specific frequency: quasicircu- 
lar initial data sequences can be used to determine the 
orbital energy at specific orbital frequencies, while full- 
NR evolutions at large separations (but for a small num- 
ber of orbits) can determine the energy flux at (nearly) 
specific frequencies. These two calculations are presum- 
ably less costly than a long numerical evolution, and they 
provide the necessary ingredients to determine the GW 
phase as a function of frequency. This scheme will be 
further explored in future work. 
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